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This  is  the  final  report  for  the  project  entitled  "Computer 
Tomography  for  Interferometric  Aerodynamic  Measurements" 
sponsored  in  the  Department  of  Mechanical  Engineering  and  Applied 
Mechanics  at  The  University  of  Michigan  by  the  Army  Research 
Office  during  the  period  20  October  1980  through  31  December 
1983  under  Contract  Number  OAAG29-81-K-0015 . 
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- 12)  Under  this  research  program  "we  have-  studied  and  developed 

several  topics  associated  with  the  use  of  multidirectional 
holographic  interferometry  with  data  analysis  by  computer 
tomography  for  application  to  complex  aerodynamic  flows.  In  this 
technique,  interferometric  data  are  recorded  virtually 
instantaneously  by  pulsed  laser  holographic  interferometry  and 
later  are  read  out  photometrically,  digitized,  and  then  processed 
by  a  computer  to  produce  density  maps  in  a  set  of  planar  cross 
sections  of  a  three-dimensional  flow.  This  program  has  included 
development  of  computer  tomography  codes  of  a  new  type, 
development  of  a  microcomputer-controlled  system  for  scanning 
interf erograms  and  preprocessing  the  data,  initial  development  of 
a  new  (non-holographic)  hybrid  optical/digital  technique  for 
recording  information  equivalent  to  that  normally  recorded  on  an 
interf erogram  of  a  flow,  the  execution  of  a  simple  prototype 
experiment,  and  considerable  progress  toward  development  of 
instrumentation  for  three-dimensional  density  measurement  in  a 
low-speed  wind  tunnel.  ^ 
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This  program  has  been  developed  on  the  premise  that 
experimental  fluid  dynamics  will  be  dominated  by  methods  which 
combine  laser  diagnostics  with  computer  control  and  analysis.  To 
date  most  work  in  the  field  has  emphasized  local,  time-resolved 
techniques  such  as  laser  doppler  velocimetry  and  various  types  of 
spectroscopy  which  have  high  accuracy  but  do  not  yield 
understanding  of  global  features  of  flows  except  by  collection  of 
vast  amounts  of  data  from  repeated  experiments.  The  current 
research  project  has  been  devoted  to  the  development  of 
techniques  which  provide  instantaneous,  whole-field  measurements 
of  aerodynamic  density,  generally  with  lower  accuracy  and  spatial 
resolution  than  the  local  techniques.  Techniques  of  these  two 
types  are  thus  complementary,  although  the  whole-field  techniques 
are  in  great  need  cf  basic  development.  Although  we  have 
considered  very  general  aspects  of  the  problem  of  measurement  of 
three-dimensional  density  fields,  the  work  we  have  conducted  is 
believed  to  have  relevance  to  efforts  such  as  that  at  the  Amy 
Aeromechanics  Laboratory  to  use  holographic  interferometry  to 
study  the  aerodynamics  of  compressible  flow  near  helicopter 
blades. 

Holographic  interf erorretry  is  the  interferometric  comparison 
of  two  optical  waves,  at  least  one  of  which  has  been  recorded 
holographically.  This  technique  can  be  used  to  produce 
inter ferograr.s  which  can  be  viewed  from  several  different 
directions.  These  interf erograms  can  be  recorded  in 
approximately  20  ns  using  a  pulsed  laser,  and  they  can  be  viewed, 


manipulated  and  analyzed  after  the  flow  experiment  has  been 
completed.  Computer  tomography  is  a  technique  for 
computationally  reconstructing  two-  or  three-dimensional  objects 
from  measured  projections,  or  line  integrals  through  them.  In 
this  case,  the  data,  are  line  integrals  of  fluid  refractive  index 
determined  by  examining  the  holographic  inter ferogram.  from 
several  different  directions.  The  result  of  such  a  computer 
tomographic  analysis  is  a  map  of  the  density  distribution  in  any 
cross  sectional  plane  of  the  flow. 

Key  topics  which  we  believed  at  the  outset  needed  research 
in  order  tc  fully  develop  the  potential  of  these  techniques  were 
augmented  during  the  research  by  other  tasks  whose  importance 
became  arparenfc  during  the  program.  These  topics  can  be 
categorized  as  Interference  fringe  Scanning  and  Processing, 
Tomographic  Algorithm  Development ,  Alternate  Technique  for 
’dave front  Recording,  and  Experimentation. 

02  122Q2QM2  RESULTS 

In  this  section  we  outline  the  key  results  obtained  duri.g 
the  research  program.  In  ost  cases,  these  results  have  been 
fully  documented  in  journal  articles,  which  are  referenced  in  th 
text  of  this  report. 

3.1  insgilerence  Fringe  Scanning  and  Processing 

Li .  1  Scanning  System 

We  have  constructed  a  system  to  record  the 
irradiance  distribution  (fringe  pattern)  of  the  waves  of  laser 
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light  reconstructed  using  holograms  of  aerodynamic  flows.  The 
system  can  be  used  to  directly  view  the  reconstructed  waves 
without  the  need  for  extra  viewing  optics.  The  image  is  recorded 
on  a  CID  array  solid  state  video  camera  in  a  i28  x  128  pixel 
format.  The  car, era  can  be  translated  through  the  image  field 
under  computer  control  by  stepping-motor  driven  precision 
translation  stages.  The  ir radiance  data  can  be  stored  on  a 
floppy  disk  and  displayed  as  an  image  on  a  CRT.  The  ir radiance 
data  from  each  pixel  can  be  directly  addressed  and  analyzed  by 
FORTRAN  programming.  An  LSI  11/23  microcomputer  is  used  for  all 
control,  recording  and  processing  functions.  The  system  permits 
one  to  calibrate  the  camera  pixel  by  pixel,  and  the  computer  will 
then  automatically  adjust  each  ir radiance  value  accordingly  when 
recording  data.  A  calibration  scheme  was  formalized,  and  a 
technical  manual  for  the  system  was  written.  The  software  and 
hardware  for  interfacing  and  controlling  the  camera,  displaying 
an  image  and  enabling  processing  of  the  data  all  were  developed 
by  graduate  students  under  the  sponsorship  of  this  ARC  contract. 


3.1.2  Exi.pg.e  Analysis 

Interferograms  of  aerodynamic  flows  are  likely  to 
have  two  undesirable  characteristics:  They  may  be  noisy,  i.e. 
contain  extraneous  irradiance  variations  and  laser  speckle,  and 
they  often  may  contain  regions  in  which  fringes  are  few  in  number 
and  therefore  broad  and  difficult  to  interpret.  Because  the 
irradiance  is  accurately  known  at  a  large  number  of  closely- 
spaced  points  in  the  field,  we  developed  a  fringe  analysis  scheme 


to  optimize  the  use  of  the  data  to  produce  accurate  phase  data 
even  in  the  presence  of  considerable  noise.  The  technique  is 
based  on  the  use  of  nonlinear  regression  analysis  to  fit 
irradiance  data  along  a  line  to  a  curve  of  the  form 

I  (x)  =  B(x)  +  E  (x)  .cos  [P  (x)  ]  , 


where  B(x)  represents  the  variation  of  background  irradiance, 

E(x)  represents  the  envelope  function  of  the  fringes  and  P(x)  is 
the  phase  variation,  i.e.  P(x)  is  the  data  needed  for  analysis  of 
the  density  field.  Although  the  approach  is  quite  general, 
polynomials  with  unknown  coefficients  were  used  represent  B(x), 
E(x)  and  P(x)  in  most  of  our  work.  The  nonlinear  regression 
algorithm  used  was  based  on  Powell's  method. 

Although  the  nonlinear  regression  analysis  requires  rather 
sophisticated  programming  and  searching,  the  approach  of  curve 
fitting  may  appear  somewhat  old  fashioned  in  comparison  with 
digital  filtering  techniques;  however,  we  used  this  approach 
because  it  recognizes  that,  physically,  the  fringe  pattern  rust 
have  an  irradiance  distribution  of  the  form  indicated  above.  We 
verified  experimentally  that  the  technique  can  be  used  to 
interpret  extren.ely  noisy  inter ferograms  with  "subfringe" 
accuracy.  The  results  of  this  investigation  have  been  presented 
in  a  journal  article  referenced  below. 
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Algorithm  Development 


We  considered  two  topics  in  the  development  of  computer 


tomographic  algorithms  for  analysis  of  holographic 
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interf erograms .  The  first  work  deals  with  the  effects  of  strong 
refraction  and  was  only  partially  sponsored  by  this  contract. 

The  second  deals  with  an  new  iterative  convolution  technique. 

1*2x1  Reconstruction  ol  Strongly  Refracting  Elides 

This  work  is  mentioned  here  briefly  because  r uch 
of  it  was  carried  out  under  different  sponsorship;  however,  it 
was  completed  under  this  project  and  has  strong  relevance  to  it. 
Virtually  all  computer  tomography  algorithms  are  based  on  the 
assumption  that  the  probing  rays  are  straight  lines  through  the 
unknown  field  so  that  the  data  represent  line  integrals,  When  a 
fluid  under  study  by  interferometry  has  strong  density  gradients, 
the  assumption  that  the  the  rays  are  straight  lines  may  not  be 
reasonable.  An  example  of  such  strong  gradients  would  be  a  shock 
wave.  We  have  developed  a  basic  algorithm  for  dealing  with  this 
case.  It  involves  an  iterative  process.  An  initial  estimate  cf 
the  field  is  made  by  reconstruction  from  the  data  without 
accounting  for  ray  bending  due  to  strong  gradients.  Path 
integrals  through  the  estimated  field  are  then  computed  by 
numerical  ray  tracing.  The  values  of  these  path  integrals  are 
then  compared  with  the  measured  experimental  data.  Based  on  the 
difference  between  the  measured  and  calculated  data,  the  estimate 
of  the  field  is  revised.  This  process  continues  iteratively 
until  an  acceptable  reconstruction  of  the  field  is  obtained. 
Details  of  this  work  are  reported  in  a  journal  article  referenced 


below. 


liie.  Xtaiative  Convolution  He t had 
Most  well-established  algorithms  for  computer 
tomography  assume  that  data  are  complete.  In  this  case 
"complete"  means  that  projections  ( interf erograms)  are  available 
over  the  entire  1?0  range  of  viewing  directions,  even  though 
they  may  be  available  only  at  discrete  intervals.  "Complete" 
also  means  that  no  opaque  objects  are  present  in  the  field  to 
block  portions  of  the  projections.  Both  of  these  conditions  are 
likely  to  be  violated  in  aerodynamic  experiments.  Physical 
constraints  ray  make  it  impossible  to  view  a.  region  such  as  the 
test  section  of  &  uird  tunnel  from  a  complete  Targe  of  viewing 
directions.  Furthermore,  an  object  such  as  the  aerodynamic  model 
under  test  i-;.  likely  to  block  a  substantial  region  of  the  field 
of  interest.  Tde  therefore*  were  particularly  interested  in 
developing  algorithms  that  would  v or k  even  when  the  data  arc* 
incomplete  ir.  either  of  the  senses  noted  above. 

The  .'.out  popular  a  Igor  i  tin  for  tomographic  reconstruction 
probably  is  the  convolution  backpro jection  algorithm  ( termed  tie 
convolution  method  here).  The  theory  of  operation  of  this 
algorithm,  is  well  understood,  and  it  runs  rapidly  and  efficiently 
on  a  computer  and  places  relatively  modest  demands  on  storage. 

The  basic  reason  for  its  simplicity  of  operation  is  that  it 
operates  in  the  object  domain.  Hnfor turately,  it  requires  that 
data  be  complete  and  available  at  equally  spaced  intervals; 
hence,  it  is  not  directly  useful  for  the  analysis  of  rn.ost 
aerodynamic  data. 

As  p:art  of  this  project  we  attempted  to  develop  an  algorithm 
that  would  retain  many  of  the  positive  attributes  of  the 
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convolution  method,  but  which  could  be  used  to  analyze  incomplete 
data  from  aerodynamic  experiments.  We  were  partially,  but  not 
completely,  successful  in  this  attempt. 

The  algorithm  we  developed  is  called  the  Iterative 
Convolution  Method.  It  operates  by  iteration  between  the  object 
and  projection  (Radon  Transform)  domains.  It  begins  by  examining 
the  projections  (fringe  data)  and  filling  in  missing  data 
according  to  some  specified  criteria.  For  example,  the 
"artificial  data"  used  to  fill  in  missing  projections  can  be 
formed  to  have  known  mathematical  characteristics  of  Radon 
transforms,  such  as  integrals  which  must  be  the  same  in  each 
projection.  Other  constraints  such  as  known  bandlimits  can  also 
be  applied  at  this  stage.  Once  this  has  been  done,  the  field  can 
be  reconstructed  by  the  efficient  convolution  method.  New 
projections  are  then  computed  as  line  integrals  across  this 
estimated  reconstruction.  Constraints  such  as  finite  spatial 
extent  can  be  applied  at  this  stage.  These  computed  projections 
are  replaced  by  measured  values  everywhere  real  data  are 
available.  A  new  reconstruction  is  then  made  by  the  convolution 
method  and  the  process  continues  iteratively. 

We  have  written  several  codes  based  on  this  genaral 
algorithm  and  have  tested  them  with  both  simulated  and  actual 
experimental  data.  Our  results  were  mixed.  When  data  are 
incomplete  because  the  entire  range  of  viewing  directions  is  not 
accessible,  the  algorithm  did  not  produce  reconstructions  as 
accurate  as  those  produced  by  other  applicable  methods  based  on 
series  expansions  of  the  unknown  field.  When  data  are  incomplete 


because  an  opaque  object  such  as  an  aerodynamic  test  model  blocks 
part  of  the  projections,  the  algorithm  appears  to  work  reasonably 
welx.  It  must  be  noted  that  this  algorithm  is  nonlinear  and  is 
being  used  to  obtain  an  approximate  solution  to  a  problem  which 
is  mathetically  ill-posed.  Its  convergence  behavior  when  the 
data  is  substantially  incomplete  is  interesting.  It  generally 
converges  toward  an  accurate  reconstruction  for  several 
iterations  and  then  diverges.  We  were  able  to  empirically 
establish  criteria  for  stopping  the  algorithm  at,  or  very  near, 
the  iteration  in  which  the  most  accurate  estimate  is  produced. 

After  initial  development  and  testing  under  this  program, 
further  work  was  done  on  this  algorithm,  including  applying  it  to 
wind  tunnel  data  by  a  post  doctoral  researcher  sponsored  by  HASA. 
Details  of  the  algorithm  and  an  analysis  of  it  are  presented  in  a 
journal  article  currently  under  review. 

3.3  £XP£XiiI£J3i£ii£I> 

We  have  been  involved  with  experiments  using  holographic 
interferometry  with  analysis  by  computer  tomography  in  three 
ways.  First,  a  simple  small  scale  experiment  was  carried  out  at 
an  early  stage  to  illustrate  the  principle.  Second,  we  have  been 
working  toward  conducting  a  low-speed  wind  tunnel  experiment  to 
measure  the  density  distribution  in  a  jet  inserted  into  a  cross 
flow.  Third,  we  have  maintained  an  active  interest  in,  and 
discussion  of,  large-scale  experiments  being  developed  at  Ames 
Research  Center.  The  first  two  of  these  will  be  described 
briefly  here. 
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3.3.1  PemonsLcatign  Experiment 

At  an  early  stage  in  this  program  we  conducted  a 
feasibility  experiment  to  show  that  the  technique  could  be  used 
to  reconstruct  a  density  field  that  was  geometrically  similar  to 
that  anticipated  in  rotating  aerodynamic  fields  such  as  the 
transonic  flows  near  helicopter  rotors.  The  actual  test  object 
was  the  temperture  field  above  a  curved  heated  wire  mounted  on  a 
rotating  table  submerged  in  water.  Data  were  recorded  by 
holographic  interferometry  at  each  of  several  rotational 
positions  of  this  object,  would  be  the  case  in  an  experiment 
in  a  rotating  system.  The  data  were  then  used  to  reconstruct  the 
temperature  field  by  using  the  convolution  method  described 
above.  The  results  were  qualitatively  correct,  and  gave  us 
confidence  to  proceed  as  planned. 

1*2 *  2  Hind  Expjiriir.gnfca.ticn 

Although  we  have  made  considerable  progress  toward 
conducting  an  in-house  wind  tunnel  experiment,  it  was  not  brought 
to  completion  during  this  contract  period.  The  problem  we  are 
working  on  is  the  insertion  of  a  high-density  gas  jet  into  an 
orthogonal  cross  flow.  The  desired  result  is  a  mapping  of  the 
density  distribution  in  crossections  of  the  jet  as  it  bends  and 
expands.  This  case  has  interest  and  relevance  to  problems  such 
as  injection  cooling  of  turbine  engines,  and  at  the  same  time  has 
many  features  in  common  with  compressible  aerodynamic  flows.  For 
example,  it  has  been  predicted  that  the  jet  will  develop  a  double 
vortex  structure  as  it  bends. 
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We  have  decided  to  create  the  high-density  jet  by  cooling  a 
gas  such  as  Carbon  dioxide  or  Nitrogen  to  a  low  temperature.  We 
have  experimented  with  various  systems  for  doing  so,  such  as  heat 
transfer  from  packed  beds  of  dry  ice  or  heat  exchange  with 
boiling  liquid  Nitrogen. 

We  have  begun  experimental  work  to  record  interf erograrr.s 
with  a  pulsed  ruby  laser,  and  have  been  working  on  a  scheme  to 
use  reflection  holography  to  record  appropriate  holograms  ir,  the 
confined  test  space  around  the  tunnel.  This  work  is  continuing 
under  a  new  ARO  contract. 

3.4  Hybrid  Processor  for  Wavefront  Recording 

Although  holography  has  many  advantages  over  other 
techniques  for  making  optical  measurements  of  aerodynamic  flows, 
it  also  has  certain  drawbacks,  primarily  based  on  the  necessity 
to  record  two  sequential  exposures.  Therefore  there  would  be 
advantages  to  a  technique  that  could  accurately  record 
quantitative  data  about  an  optical  wavefront  that  has  been 
distorted  by  passage  through  a  flow  without  the  need  for  two 
recordings.  We  therefore  explored  the  use  of  hybrid  processing, 
i.e.  the  combining  of  optical  and  computer  processing  of  optical 
data. 

The  system  we  developed  is  based  on  the  use  of  a  simple, 
well-known  technique  of  optical  Fourier  processing.  The  wave 
front  to  be  studied  enters  a  processor  that  consists  of  two 
lenses  seperated  axially  by  the  sum  of  their  focal  lengths.  At 
the  back  focal  plane  of  the  first  lens  the  resulting  irradiance 
pattern  is  known  to  represent  the  Fourier  transform  of  the 
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entering  wave.  In  this  transform  plane  we  place  a  filter  that 
differentiates  the  signal,  i.e.  the  output  of  the  second  lens  is 
the  derivative  of  the  input  wave.  This  system,  is  very  similar  to 
a  classical  schlieren  system  with  laser  illumination,  except  the 
filter  is  more  complicated  than  a  knife  edge. 

The  output  of  the  optical  processor  is  recorded  by  a  solid 
state  array  car. era  and  stored  in  a  computer  where  it  is 
integrated  to  give  the  phase  distribution  across  the  field.  This 
phase  distribution  is  the  information  we  normally  would  record  by 
interferometry.  The  processing  described  so  far  is 
straightforward  and  well  known.  However,  such  coherent  optical 
processing  is  notoriously  noisy,  so  much  so  that  it  is  not  very 
useful  for  practical  applications  to  flow  visualization  and 
measurement.  We  have  developed  a  technique  for  greatly  improving 
it.  In  our  technique,  the  computer  records  three  different 
outputs  from  the  optical  processor.  Fach  corresponds  in  essence 
to  a  slightly  different  position  of  the  filter  in  the  Fourier 
plane.  When  these  three  signals  are  combined  in  a  particular 
algebraic  manner  in  the  computer,  most  of  the  noise,  including 
that  caused  by  nonuniform  background  illumination  ar.d  much  of  the 
speckle  noise  are  eliminated.  We  believe  that  this  hybrid 
processing  technique  is  a  potentially  important  flow  diagnostic 
tool.  This  work  is  described  in  detail  in  a  journal  article 
referenced  below. 
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Hybrid  processing  for  phase  measurement  in  metrology  and 
flow  diagnostics 


Ivan  Prikryl  and  Charles  M.  Vest 


A  hybrid  <>ptic.il  dt.iMl  pr<Ke.*.*ing  -thrme  u>r  measuring  phase  distrilnitionj*  i>  described  and  drm<<n 
pirated.  It  :*  tntended  t»*  lie  an  alternative  t«»  intertert»meinc  methods  of  measuring  uptical  path  length 
changes  in  ?l«*w  tii.iu'n"*tu*  and  ian  also  Kt*  used  as  a  tlow  visualization  technique  The  processing  *»  hetne 
enables  i»ne  to  make  accurate  measurements  **i  phase  at  arbitrary  points  in  the  image  plane  The  *wt*m 
is  based  «»n  i  simple  coherent  optical  Fourier  processor  but  incorporates  three  separate  measurement*  m*i 
postdetection  digital  processing  to  eliminate  extraneous  parts  o|  the  *ignal.  The  addition  oi  a  hoio.raphu 
filter  to  the  ^\*tem  enables  one  t*»  measure  detormation  or  displacement  oi  dittuselv  reileiting  opaque 
objects.  The  'eihmque  is  demonstrated  by  using  it  to  vi>udliie  the  flow  of  an  expanding  compressible  ..a* 
jet  and  to  measure  the  optical  path  length  through  a  heated  plume  ot  air 


I.  Introduction 

Many  physical  measurement  techniques  require  the 
determination  ot  the  spatial  distribution  of  phase  of 
coherent  light  in  some  plane.  Because  detectors  are  not 
capable  of  responding  at  optical  frequencies,  phase 
distributions  must  be  coded  into  irradiance  distribu¬ 
tions  which  can  be  detected  or  recorded  by  square- law 
detectors  such  as  photosensors  or  photographic  film. 
The  most  common  technique  for  such  encoding  is  in¬ 
terferometry,  whereby  the  phase  distribution  is  dis¬ 
played  as  a  fringe  pattern.  F rom  such  a  pattern,  rela- 
tive  phases  at  the  center  of  fringes  can  be  determined. 
However,  when  accurate  measurements  must  be  made 
at  other  points,  such  as  on  a  uniform  grid  of  points 
which  do  not  coincide  with  fringe  centers,  additional 
effort  is  required.  For  example,  a  heterodyne  system 
can  be  used.1  or  irradiance  measurements  at  the  points 
can  be  determined  by  interpolation.-1  This  interpola¬ 
tion  is  not  straightforward  because  the  desired  phase 
distribution  is  the  argument  of  a  fringe  function,  most 
often  sinusoidal. 

In  automated  systems,  irradiance  distributions  can 
be  recorded  directly  using  vidicons  or  solid-state  de¬ 
tector  arrays,  thus  bypassing  the  intermediate  step  of 
photographing  fringe  patterns.  In  this  paper  we  point 
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out  that,  when  the  capability  of  accurately  and  rapidly 
recording  irradiance  distributions  exists,  it  may  be  de¬ 
sirable  to  reconsider  other  ways  of  mapping  phase  dis¬ 
tributions  into  irradiance  variations,  namely,  coherent 
optical  Fourier  processing  and  related  filtered  schiieren 
and  shearing  interferometric  techniques.  In  particular, 
we  discuss  ways  of  recording  phase  gradients,  w  hich  can 
be  integrated,  if  necessary,  to  determine  phase  distri¬ 
butions.  We  first  describe  the  approach  in  a  form 
suitable  for  applications  in  which  a  plane  object  wave 
has  been  slightly  distorted  by  passage  through  a  tran>- 
parent  medium  or  reflected  by  a  specular  surface.  This 
would  be  applicable  to  wind  tunnel  diagnostics,  heat 
transfer  studies,  plasma  diagnostics,  etc.  An  extension 
of  the  technique  will  then  be  described  in  which  a  ho¬ 
lographic  filter  can  be  used  to  measure  phase  differences 
due  to  distortion  or  mo'  on  of  diffusely  scattering  sur¬ 
faces. 

II.  Method 

In  the  proposed  method  the  phase  distribution  of  a 
slightly  distorted  plane  wave  of  coherent  light  is  mea¬ 
sured  by  passing  the  wave  through  a  simple  optical 
processor  which  performs  a  total  differentiation.  The 
output  of  the  processor  is  recorded  with  a  vidicon  or 
solid-state  array  camera  and  subsequently  analyzed 
further  by  a  digital  computer.  It  will  be  demonstrated 
that  this  hvbrid  approach  alleviates  many  of  the  prac¬ 
tical  difficulties  associated  with  coherent  optical  pro¬ 
cessing.  Such  a  system  is  show  n  in  Fig.  1. 

Let 

u  •*  t  .“.i  *  4'U  »  expjo  *».t  .\  *| 

be  the  complex  amplitude  in  a  plane  :  *  constant  of 


»UM 


Fi*.  1.  Optical  pr»H.essor.  The  tiller  is  placed  in  the  back  h>cal  plane 
i>l  the  transtorminK  lens.  In  general,  the  output  lens  is  aligned  to  form 
an  image  of  the  test  region  in  the  detector  plane. 


interest.  Our  objective  is  to  determine  oix.x).  As 
shown  in  Fig.  1,  the  processor  consists  of  a  Fourier 
transforming  lens,  a  Fourier  filter  in  its  back  focal  plane, 
and  a  second  lens  which  reforms  the  processed  wave 
front  so  that  its  irradiance  distribution  can  be  recorded 
photographically  or  electronically.  The  Fourier  filter 
used  in  the  current  study  is  designed  so  that  the  pro¬ 
cessor  forms  the  total  differential  of  the  input  signal. 

Now  suppose  that  the  filter  that  performs  the  dif¬ 
ferential  operation  is  translated  laterally  (parallel  to  the 
X  axis)  a  distance  - tf/k ,  where  f  is  the  focal  length  of 
the  transforming  lens  and  k  =  2jt/A  is  the  wave  number. 
The  translated  filter  will  differentiate  the  function 

ut.Y.Y)  =  vi.Y.Y'l  expli (^i A'.  Vl  +  <A’|l  i2) 

instead  of  the  function  U)  and  simultaneously  will 
multiply  the  differential  by  the  value  exp(  —itX ).  The 
output  of  the  processor  will  be 


nonzero.  We  note  that  the  use  of  three  measurements 
to  remove  ambiguity  and  to  remove  extraneous  com¬ 
ponents  of  the  irradiance  distribution  is  conceptually 
similar  to  certain  techniques  of  two-reference  beam 
holographic  interferometry,*  real-time  holographic 
interferometry,4  and  interferometric  real-time  phase 
measurement.* 

Having  determined  the  total  differential  d&iX.Y) 
over  the  X-Y  plane  by  using  Eq.  (7),  the  phase  d>iP)  at 
any  point  P  can  be  found  by  path  integration  along  any 
curve  which  connects  P  and  any  point  Pti  where  the 
phase  has  a  known  value  0(Pn): 

o\P)  =  doiP  )  +  i*i 

Such  integration  of  optically  differentiated  phase  dis¬ 
tributions  has  been  reported  by  Sprague  and  Thomp¬ 
son*  who  used  an  optical  integration  technique. 

If  the  wave  front  under  consideration  contains  an 
aberration  &niX,Y)  whose  order  does  not  exceed  that 
of  the  variation  0LY,V)  under  study  the  technique  can 
still  be  used.  In  this  case  the  complex  amplitude  Eq. 
(1)  contains  the  0(X,YP)  +  OnlX.V)  as  its  phase  term, 
where  onlX.Yr)  may  be  due  to  misalignment  .mperfect 
optical  elements,  or  low  quality  test  section  windows. 
If  0niX.Y)  does  not  change  with  time,  two  sets  of  three 
exposures  can  be  recorded — one  set  before  the  change 
o(X.Y')  is  introduced  and  one  set  after  it  is.  Because 
the  computed  phases  are  stored  in  a  computer,  they  can 
be  subtracted  to  give  information  only  about  0(X,Y'). 
This  is  illustrated  experimentally  in  Sec.  IY\ 


d  m.Y.Yi  =  Scki.Y.Vt  +  t|H.Y.Vi|<foi.Y.Vi  +  (dX\\ 
•  exp|fO<.Y.Y'>|. 
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The  irradiance  distribution  of  the  output  would  be 


/i.Y.Vi  *  (c/iM.Y.l  )|-  +  g-i.Y.Vi|rfoi.Y.yi  +  <dX\- 

This  irradiance  distribution  is  recorded  electronically 
and  stored  in  the  memory  of  a  digital  computer  for 
subsequent  additional  processing. 

The  total  differential  d<MX,Y')  which  we  desire  to 
measure  is  present  in  this  expression  for  the  irradiance 
HX.Y)  but  must  be  separated  from  the  effects  of 
g'-lX.Y)  and  [d£(X,  Y')J-  Such  a  separation  can  be 
made  if  three  different  irradiance  measurements  are 
made,  each  corresponding  to  a  different  position  of  the 
filter.  In  particular,  let  the  measurements  /0(X,Vh 
/ i(X.Y'),  and  I>  (X.V)  correspond  to  values  *  3  0,  (  * 
-fo*  and  f  =  fo,  respectively.  From  Eq.  <6)  it  is  easily 
shown  that 


doi.Y.Vi  *  -  Unrf.Y > 


htX.Yi  -  / [ i A'  >  ) 

/ , i A' .V»  +  /ji.V.Vi  -  j/..t A\y»  * 


f  7  > 


Note  that  the  denominator  of  Eq.  (7)  equals  m2g-\X,Y) 
<MdX)~  which  is  nonzero  as  long  as  g(X.Y')  and  f,»are 


III.  Differential  Filters 

We  mention  here  two  types  of  filter  which  could  be 
used  in  the  optical  processor  described  in  the  previous 
section.  The  first  is  the  ideal  Fourier  filter  for  forming 
a  differential,  and  the  second  is  a  composite  grating  that 
approximates  the  desired  filter. 

The  ideal  filter  to  perform  a  differentiation  is  based 
on  the  differentiation  theorem  of  Fourier  transforms. 

i  a.Y ;  =  :uF*u.i  ».  1 
*)Af\  A  .  V1  ■  u  F*u  .i  1 


where  F(u.r)  =  J|/(X,Vl|  is  the  2-D  Fourier  transform 
of  fiX,Y).  For  our  purposes,  fiX.Y)  represents  the 
complex  amplitude  of  the  input  field  to  the  processor. 
The  first  lens  of  the  processor  produces  the  Fourier 
transform  of  this  field  in  its  back  focal  plane.  If  £  and 
7  are  the  physical  coordinates  in  this  plane,  the  complex 
amplitude  there  is 


/’i. <*i)//i|.  ii«u 

where  k  *  2sr/A  and  {  is  the  focal  length  of  the  trans¬ 
forming  lens. 

According  to  Eqs.  (9)  the  filter  must  have  an  ampli¬ 
tude  transmittane  which  varies  linearly  in  the  transform 
plane  and  must  produce  a  phase  shift  to  simulate  neg¬ 
ative  values  of  u  and  c.  Hence  the  filter  must  have  the 
following  intensity  transmittance: 

*  C-ik*  7  +  ill* 
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and  havea  \/2  thin  film  layer  in  the  region  ik/f)  (£  + 

<  0  in  order  to  introduce  the  phase  shift 

(l2, 

1  <T  lor  \k  f  +  rp  <  0. 

The  constant  C  should  have  a  value  such  that  the  in¬ 
tensity  transmittance  is  unity  at  the  outer  edge  ot  the 
filter: 


C-ifc£  /  +•  W  -  1  1  l-?l 

The  processor  using  this  filter  will  produce  a  single 
output  wave  of  the  form  of  Eq.  to)  where  the  parameters 
in  Eqs.  l3)-i5)  must  be  interpreted  as  follows: 

i/io  X.  V»  •-(  .Y . V  *  ,.Y  > -v. * .Y 
it. Y 

■iV  •  -•(' 

Then  Eqs.  (6)  and  (7)  are  also  valid  if  we  take 

uiAV  -  •(/>'»-  =  -  <’  * 

One  possible  realization  of  the  ideal  filter  for  1-D 
differentiation  has  been  used  by  Cody-  to  visualize 
gradients  in  compressible  tlows.  In  his  experiments  two 
separate  filters  were  used — one  with  a  quadratic  in¬ 
tensity  transmittance  and  one  with  a  step  change  of 
phase.  These  filters  were  placed  at  conjugate  planes 
of  an  additional  lens  in  the  processor. 

Other  ways  of  forming  differentials  without  using  a 
\/2  thin  film  have  been  investigated.  Usually  some 
compensation  must  be  used  to  correct  for  nonlinear 
relationships  between  the  detected  irradiance  and  the 
derivative  of  the  phase.  This  compensation  may  in¬ 
volve  heavy  postdetection  electronic  processing,'  strong 
predetection  biasing,*  or  the  use  of  a  square  root  filter 
rather  than  a  linear  filter.1'1 

A  composite  grating  can  he  used  as  a  Fourier  filter  to 
perform  the  operation  in  Eq.  (3)  on  a  more  approximate 
basis.10  An  appropriate  composite  grating  can  be 
produced  very  easily,  and  one  was  used  in  the  illustra¬ 
tive  experiments  described  in  Sec.  IV. 

We  first  consider  a  composite  grating  filter  which 
performs  an  approximate  differentiation  in  one  di¬ 
mension.  It  consists  of  two  superimposed  sinusoidal 
gratings  which  have  slightly  different  frequencies  and 
which  are  shifted  by  half  a  period  relative  to  each  other. 
Each  grating  will  produce  three  component  waves  in  the 
Fourier  plane — the  central  order  and  two  diffracted 
orders.  The  propagation  angle  of  the  diffracted  waves 
is  proportional  to  the  grating  frequency.  Hence  each 
diffracted  order  leaving  the  composite  grating  will 
contain  a  pair  of  waves  having  slightly  different  prop¬ 
agation  directions  and  a  relative  phase  shift  of  tt. 
Correspondingly,  in  the  output  I  inverse  Fourier  trans¬ 
form)  plane  each  order  will  consist  of  two  slightly 
sheared  and  shifted  waves.  Hence  the  processor  with 
this  filter  acts  as  a  modified  shearing  interferometer.  If 
T  is  the  period  of  one  grating  and  the  period  of  the  other 
grating  is  T  +  AT,  the  shear  will  be 

A.Y  -  \n±Ti  r-.  -ini 


where  f  is  the  focal  length  of  the  first  lens  of  the  pro¬ 
cessor. 

Stated  differently,  the  impulse  response  of  the  pro¬ 
cessor  is 


Ai A’  +  AAWi  -tiiA.Vi,  •  i“i 

where  h  is  the  Dirac  delta  function,  i  In  this  discussion 
we  assume  that  the  origin  of  the  coordinate  system  has 
heen  shifted  to  the  center  of  one  of  the  diffracted  or¬ 
ders.)  So.  if  ui  .V,  V)  is  the  input  signal,  the  output  is 

in  A  +  A  A . Vi  -  u« X. Yt  <  I  *o 

If  AX  is  sufficiently  small  to  express  the  difference  in 
Eq.  1 18)  using  only  the  first  two  terms  of  a  Taylor  >eries. 
we  see  that  the  output  is  approximately 

•*«  YY  i  A  A  ;  m 

If  the  wave  having  complex  amplitude  utX.V)  given 
by  Eq.  (2)  is  acted  on  by  this  processor,  the  output  will 
he 
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The  corresponding  irradiance  distribution  is 


/. A.Vi  =  if-i.Y.Vi  «-  uiA.V*  <►  —  AX 

l  '.Y  I 

-iM.V.Vi  L'I.Y.Vi  +  -i.Y 

l  .‘.v 
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Finally,  if  the  operation  described  in  Sec.  II  were  carried 
out.  the  expression  corresponding  to  Eq.  (7)  would  be 


1  -  r.»»e,-AA)  /_■«  A.Vi  -  /;<A.V» 
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Hence  if  AX  and  f,,AX  are  sufficiently  small,  then 

A  Vi  I  t  /  .  A.  V'  -  /  •  A.  V» 

- A  A  -  -  *#  A.Y  - - - -  . 

•  *A  J  /  » A.> »  W_|.\.Y>  -.7  .  \  >  ■ 

which  is  of  the  desired  form  un  one  dimension!.  In 
particular.  AX  must  he  >mall  enough  so  that 

UMi’.i-  ■•A*AA|  '  fi.Y  'A.Y  ■  ;u 

throughout  the  field  for  Eq.  ( 23)  to  he  valid. 

As  would  he  expected  on  physical  grounds.  Eq.  *22i 
will  he  unambiguous  only  if 


is  a  criterion  which  must  he  met  by  the  processor  in 
order  to  use  Eq.  1 22).  Analysis  also  indicates  that  ;hr 
angle  in  Eq.  r22)  must  he  chosen  such  that 
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*  t  <  *  . A  <  t .  ( J7> 

Preferably,  the  value  of  f,,AA'  should  be  chosen  to  be 
:  2.  If  simultaneously  we  choose  a  value  of  AX  close 
to  the  limit  in  Eq.  (26),  we  obtain  the  maximum  possible 
difference  between  the  intensities  /_>  and  / 1.  This  dif¬ 
ference  varies  from  zero  i  for  do  =  0)  to  the  maximum 
intensity  which  can  be  reached  by  constructive  inter¬ 
ference  in  the  instrument.  The  best  value  for  (a  is 

.AL,,.  •>» 

To  form  an  approximate  total  differential  as  required 
in  Eq.  (7i,  the  composite  grating  can  consist  of  three 
sinusoidal  amplitude  gratings  which  yield  the  impulse 
response 11  ’ 

AV  *  AAW*  -  '*  \\>  ■ 

This  filter  will  be  subject  to  restrictions  analogous  to  Eq. 
(24b  Also,  the  radius  R  of  the  filter,  which  should  also 
he  the  aperture  stop  of  the  processor,  is  such  that 

-  f<  t  »  \  _,A.\  «■  -•  . 

where  /  is  the  focal  length  of  the  transforming  lens. 
Condition  (20)  permits  the  impulse  response  to  be  ap¬ 
proximated  by  delta  functions  as  in  Eq.  «29>. 

Note  that  if  composite  grating  filters  are  used  they 
must  be  designed  so  t  hat  in  a  given  experiment  AX  and 
AV  will  meet  the  criteria  which  depend  on  the  maxi¬ 
mum  value  of  the  phase  distribution  ■;*» X.  V». 

IV.  Experimental  Results 

Two  experiments  were  conducted  in  order  to  dem¬ 
onstrate  this  technique,  fn  the  first,  the  optical  pro¬ 
cessor  was  used  to  visualize  a  compressible  gas  jet.  In 
the  second,  the  hybrid  processor  was  used  to  measure 
the  temperature-dependent  optical  path  length  distri¬ 
bution  across  a  natural  convective  plume.  In  the  latter 
case  the  ability  to  remov  e  the  effects  of  wave  front  ab¬ 
errations  was  also  demonstrated. 

Both  experiments  were  carried  out  using  a  composite 
grating  filter  constructed  by  recording  interference 
patterns  on  a  holographic  plate.  The  grating  performed 
an  approximate  1-D  differential  operation.  It  produced 
a  background  wave  and  two  composite  diffracted  waves 
traveling  at  angles  of  ±2. >°  w  ith  respect  to  the  g  axis  and 
having  linear  shears  AX  =  0.2  mm  in  object  space.  The 
impulse  response  of  either  off-axis  component  is  given 
by  Eq. (17b 

In  the  first  experiment  the  object  was  a  jet  of  com¬ 
pressed  nitrogen  gas  flowing  from  a  nozzle.  In  this  ca^e 
the  output  of  the  optical  processor  was  recorded  pho¬ 
tographically.  Figure  2  is  a  photograph  of  the  nozzle 
with  no  flow.  Note  the  bright  lines  which  form  because 
of  the  steep  horizontal  gradients  in  the  image  at  the 
vertical  edges  of  the  nozzle.  The  width  of  these  lines 
is  defined  by  the  shear  AX  *  0.2  mm.  Figure  2  is  a 
photograph  recorded  while  the  gas  jet  was  present.  The 
edges  of  the  jet  and  its  periodic  structure  are  clearly 
visible. 

In  the  second  experiment  the  output  of  the  optical 
processor  was  recorded  by  using  a  solid-state  (TD  array 


Fig.  Processed  image  nt  a  nozzle  with  nu  il<-u 


Fig.  Professed  image  * j  •«  •  mtfogen  gj>  i)i tn <m  a  n> *zz)e 
A  composite  grating  tilter  is  u^ed  to  ditterentute  approximatelv  m 
’he  horizontal  direction 

camera  with  128  X  128  pixels.  The  images  recorded  bv 
the  camera  were  stored  in  and  processed  bv  a  micro¬ 
computer.  The  processing  consisted  of  two  parts. 
First,  each  signal  was  averaged  over  a  few  pixels  in  the 
vertical  direction  to  suppress  noise  introduced  by  the 
photographic  grating.  Second,  after  the  three  required 
images  were  recorded.  Eq.  (22)  was  evaluated.  Finally. 
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the  result in^  values  of  the  path  length  gradients  were 
integrated  to  obtain  valuer  of  the  path  length  itself. 
Two  sets  of  measurements  were  made  in  this  way  and 
were  subtracted  to  give  the  desired  phase  distribu¬ 
tion. 

The  light  wave  was  tilted  slightly  in  order  to  introduce 
an  extraneous  phase  -hift  which  was  considered  to  be 
noise  in  the  signal.  Figure  4  is  a  plot  of  the  optical  path 
length  distribution  measured  in  this  way.  The  dashed 
curve  is  the  optical  path  legth  distribution  due  to  the 
-urn  of  the  noise  and  the  heated  plume.  The  solid  curve 
is  the  distribution  after  the  noise  was  filtered  out  by 
subtracting  the  measurement  made  without  the  heated 
plume  from  that  made  with  the  plume. 

V.  Application  to  Complicated  or  Diffuso  Waves 

In  Secs.  II— IV  we  have  considered  the  application  of 
a  hybrid  processor  to  measure  the  phase  distribution  of 
a  relatively  simple  wave  front,  e.g..  a  plane  wave  which 
has  traveled  through  a  wind  tunnel  test  section  or  has 
been  distorted  slightly  by  a  reflective  or  refractive  op¬ 
tical  element.  If  the  wave  of  interest  is  diffused  prior 
to  passing  through  a  transparent  object  or  is  scattered 
by  a  diffuse  opaque  object  whose  deformation  is  to  be 
measured,  the  system  can  be  used  if  a  holographic  filter 
is  added  to  it. 


The  required  holographic  system  is  shown  in  Fig.  5. 
The  diffuse  object  surface  is  imaged  onto  a  recording 
medium.  In  the  recording  plane  z‘  =  0,  the  complex 
amplitude  of  the  object  wave  is 

u..u  .>  \  *  ao  exp(i4»i  r  .>  >|.  '-U» 

An  off-axis  image  plane  hologram  is  recorded  by  adding 
to  Uni  x  .y  )  a  reference  wave 

UK u  »  *  or  exp[Ma.r  '  +  v  I.  1  -iji 

This  hologram  is  developed  in  place  or  is  developed  and 
precisely  replaced  in  the  system  after  development.  It 
has  an  amplitude  transmittance  of  the  form 

t ...  *  t  -iv  «  +  .i<rxp|rfe».v  •  -  n.r  ' 

-  r-\p:-  Kv  -  :)  ■  '  '* 

The  disturbance  of  interest  is  then  introduced  so  that 
the  object  wave  at  the  hologram  plane  becomes 

u  -i  .  •  *  e\pt;  i  *!>•  r  v  >11.  4- 

where  the  distribution  of  phase  change  ’»  is  the 

quantity  to  be  mea>ured. 

If  u.  \x  ,\  ’)  alone  illuminates  the  hologram  described 
by  Eq.  «.W>.  the  output  of  the  hologram  will  be 

u  •  \  r  i  . ■  '  =  expi i'  +  ■:■»! 

♦  i  e\pi<‘ '  -  -  •  t  ‘I 

+  J  e\pi*  «■.!+•  i«.\  •  *  > 

Equation  indicates  that  the  wave  traveling  away 
:r< >m  the  hologram  in  the  same  direction  as  the  original 
reference  wave  has  the  desired  phase  structure  and  can 
-erve  as  the  input  wave  for  the  hybrid  processor  < see  Fig. 
•>i.  The  complicated  structure  h  which 

remains  constant  between  exposures,  has  been  removed 
from  this  wave.  Note  that,  if  the  :  axis  of  the  coordi¬ 
nate  -v-tem  used  in  writing  Eq.  1 1  >  is  taken  to  be  in  the 
propagation  direction  of  the  holographic  reference  wave 
the  last  term  in  Eq.  i is  identical  to  that  in  Eq. 
<  l »  with 


If  the  motion  < >r  phase  disturbance  creates  a  s\  'tern 
of  fringes  localised  otf  the  surface  of  the  object  under 
-tudy.  the  imaging  -v-tem  should  be  focu>ed  **n  *he 
region  of  loeaii/.at;*  A  «-r  the  aperture  must  hf-ci- 
ticiently  small  to  torn;  a  w-  ’-defined  structure  ot  pha-e 
difference  .'i.v.y  »  m  the  hologram  plane.  Bv  placing  the 
aperture  stop  in  the  focal  plane  of  the  imaging  lens  it  can 
he  insured  that  the  surface  of  localization  is  on  the 
object  surface.  : 

If  necessary,  the  wave  leaving  the  holographic  f  ilter 
can  be  recorded  at  any  instant  of  time  on  a  second  off- 
axis  hologram.  This  wave  can  be  reconstructed  at  any 
later  time  and  processed. 

VI.  Conclusion 

In  summary,  we  have  described  and  demonstrated 
techniques  for  using  hybrid  processing  to  visualize  and 
measure  phase  changes  due  to  variations  of  refractive 
index  in  transparent  media  or  deformation  or  dis¬ 
placement  of  opaque  objects.  Potential  advantages  of 


2848  APPLIED  OPTICS  Vol.  22.  No.  10  15  September  1983 


this  technique  relative  to  interferometry  include  the 
ability  to  make  accurate  phase  measurements  at  arbi¬ 
trary  points  (not  just  at  fringe  centers)  and  the  ability 
to  make  such  measurements  without  ambiguity  in  sign 
of  phase  change.  The  technique  uses  digital  processing 
to  remove  much  of  the  optical  noise  associated  with 
coherent  optical  processing. 


The  authors  wish  to  acknowledge  the  assistance  of 
John  B.  Schemm  who  developed  the  camera/computer 
system  used  to  perform  the  digital  processing.  This 
work  was  sponsored  by  the  L\S.  Army  Research  Of¬ 
fice. 
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Fringe  pattern  recognition  and  interpolation  using 
nonlinear  regression  analysis 


John  B.  Schemm  and  Charles  M.  Vest 


l, east -square  error  criteria  are  used  to  lit  l-D  interference  fringe  pattern  irradiance  data  to  a  physically 
meaningful  function  of  the  form  Hx)  =  Hix\  +  Eix)  cos|Plxl|.  where  Bix).  E  lx),  and  Pit)  are  low -order 
polynomials.  I  his  procedure  is  intended  to  complement  digital  fringe  recognition  by  providing  a  method 
tor  smoothing  and  interpolating  among  fringe  position  data  when  the  number  of  Iringes  is  small,  there  ire 
more  than  ten  irradiance  measurements  per  fringe,  and  accurate  phase  values  are  needed  at  arbitrary  lina- 
tions  in  the  field. 


I.  Introduction 

Accurate  digital  recognition  and  interpolation  of  in¬ 
terference  fringe  patterns  in  the  presence  of  laser 
speckle  and  extraneous  slowly  varying  background  ir¬ 
radiance  are  important  for  development  of  automated 
systems  for  interpreting  interferograms.  The  infor¬ 
mation  of  interest  for  quantitative  analysis  of  an  in- 
terferogram  is  the  phase  distribution  encoded  in  the 
fringes  or  a  derivative  of  this  phase  at  arbitrary  positions 
within  the  fringe  pattern.  If  fringe  spacing  is  large,  it 
is  necessary  to  interpolate  among  fringe  position  data 
to  find  this  phase  or  its  derivatives  at  such  arbitrary 
positions.  It  is  frequently  necessary  to  smooth  the 
fringe  position  data  (particularly  if  there  is  noise  such 
as  speckle)  to  achieve  accurate  results.  In  this  paper 
we  present  a  method  for  performing  this  smoothing  and 
interpolation  that  is  intended  to  complement  digital 
fringe  recognition. 

Fringe  recognition  generally  involves  locating  the 
centers  of  dark  and  light  fringes.  Recently,  various 
researchers  have  developed  schemes  to  automate  this 
process.  Video  or  solid  state  array  cameras  provide 
irradiance  measurements  on  a  matrix  of  points  (pixels) 
over  the  interferogram.  Image  enhancement,  level 
slicing,  edge  detection,  and  related  processing  tech¬ 
niques  are  then  used  to  locate  fringe  centers  and  to  track 
them  along  continuous  fringes  (see  e.g..  Ref.  1). 

The  authors  and  their  colleagues  frequently  work 
with  holographic  interferometry  of  transparent  media. 
In  this  work,  and  in  many  other  applications  of  inter¬ 
ferometry.  the  number  of  fringes  in  a  pattern  may  be 
very  small.  In  such  cases  there  is  frequently  not  enough 
information  for  accurate  evaluation  of  the  phase  dis¬ 
tribution.  There  are  not  enough  fringe  centers,  and 
their  positions  are  not  known  with  sufficient  accuracy. 
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However,  each  individual  irradiance  measurement 
carries  information  about  the  phase  of  the  fringe  at  its 
position.  The  method  reported  in  this  paper  uses  all 
the  information  available  in  the  irradiance  measure¬ 
ments  to  smooth  and  interpolate  the  phase  of  the  fringes 
accurately.  This  is  accomplished  by  using  nonlinear 
regression  analysis  to  fit  an  algebraic  expression  to  the 
full  irradiance  field. 

To  develop  this  technique,  l-D  fringe  patterns  were 
digitized  as  sets  of  values  of  irradiance  measured  at 
equally  spaced  intervals.  Least-square  error  criteria 
were  used  to  fit  these  measured  values  to  a  function  of 
the  form 

h  x  i  =  B<  x  i  +  Em  cos  (Pi  x  i|.  ( 1  i 

where  Bix).  E\x).  and  Pix )  are  low-order  polynomials. 
fl(.v)  describes  the  variation  of  the  background  irra¬ 
diance  due  to  such  things  as  nonuniform  object  illumi¬ 
nation.  Eix)  is  an  envelope  function  related  to  the 
visibility  of  the  fringes.  Pl.t)  is  the  phase  of  the  fringes 
and  is  usually  the  quantity  of  interest. 

This  method  is  intended  to  circumvent  some  of  the 
problems  encountered  by  most  fringe  recognition 
schemes  due  to  the  undesirable  characteristics  of  irra¬ 
diance  measurements.  Two  particular  problems  are 
noise  and  nonuniform  illumination.  Noise  may  be  due 
to  laser  speckle  or  thermal  and  electrical  effects  in  the 
sensor  system.  Fringe  brightness  and  visibility  may 
vary  across  the  image  because  of  fringe  localization, 
nonuniform  object  illumination  or  reflectivity,  a  non- 
uniform  holographic  reference  wave,  a  nonuniform 
holographic  reconstruction  wave,  or  any  combination 
of  them. 

II.  Discussion 

Automated  fringe  pattern  recognition  systems  have 
tended  to  copy  the  operation  of  the  eye.  Such  systems 
recognize  fringes  and  represent  them  in  binarv  form  as 
bright  or  dark  regions.  Variations  in  brightness  and 
visibility  generally  are  ignored.  This  makes  accurate 
estimates  of  the  phase  distribution  across  fringes  dit- 
ficult. 
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When  two  monochromatic  linearly  polarized  wave 
fronts  with  complex  amplitudes  U \  and  IK 

UiU>  3  .4{U>  exp[-H*>.u>|.  »-* 

L’jit*  *  exp(— io-jCt  *|.  <;0 

interfere,  the  irradiance  Iix )  along  a  line  through  the 
field 

/i.t»  3  j U t x » | -  *  1Ui»x)  +  U v * | -  Ui 

has  the  form  of  Eq.  ( 1).  When  one  is  describing  the 
form  of  interference  fringes  it  is  common  to  assume  that 
fit. rb  the  background  irradiance.  and  Eix  >.  the  envelope 
function,  are  constant  and  equal.  This  is  not  a  rea¬ 
sonable  assumption  for  automated  digital  processing 
of  interference  fringes. 

Computerized  video  systems  measure  irradiance  /, 
at  a  large  number  of  positions  The  irradiance 
measurements  /,  have  optical  and  electrical  noise  that 
make  it  difficult  to  evaluate  the  phase  of  interference 
fringes  directly  from  them.  In  our  approach  an  alge¬ 
braic  expression  l(xt)  is  fitted  to  the  set  of  irradiance 
measurements  by  least -square  error  methods.  A 
physically  meaningful  form  for  lix,  >. 

/it.  i  *  fU  r.  i  £<  r  »  cosj/h  v  »|.  i*) 

has  been  used  with  low-order  polynomials  for  Bi. r  K 
),  and  Pixt ).  If  e,  is  the  difference  between  an  ir¬ 
radiance  measurement  and  its  fitted  value, 

e,  =  / 1  x. )  -  / IT.  ».  oU 

and  the  sum  of  the  squares  of  the  errors  is  denoted  bv 

E . 

E  *  ^  t7* 

the  coefficients  in  the  estimate  /(*,)  are  selected  to 
minimize  E. 

In  this  manner  the  phase  function  P{x)  is  approxi¬ 
mated  by  the  function  Pixt )  with  the  constraint  that 
l(xt)  describes  the  modulated  cosine  fringe  function 
expected  on  physical  grounds.  The  phase  function  is 
then  available  for  interpolation,  calculating  derivatives, 
or  any  other  type  of  evaluation. 

III.  Difficult)** 

The  two  major  drawbacks  of  this  method  are  both 
related  to  use  of  nonlinear  regression  analysis  and 
suggest  that  it  must  be  applied  with  great  care.  All 
nonlinear  regression  algorithms  are  iterative.  As  such 
they  require  significant  amounts  of  computer  time. 
Also,  unlike  linear  regression  analysis,  they  require 
initial  guesses  of  the  values  of  all  coefficients. 

Nonlinear  regression  analysis  algorithms  converge 
to  the  local  minima  of  E  nearest  the  initial  guess.  This 
may  not  be  the  correct  solution.  In  the  present  appli¬ 
cation  the  irradiance  estimate  is  sufficiently  nonlinear 
that  there  can  be  a  large  number  of  these  local  minima. 
For  the  procedure  to  converge  to  the  correct  description 
of  the  fringe  pattern  the  algorithm  must  have  very  good 
quality  initial  guesses. 


Good  initial  guesses  can  he  made  by  using  a  fringe 
recognition  scheme  to  locate  approximately  the  position 
and  magnitude  of  the  fringe  maxima  and  minima.  The 
fringe-order  number  can  then  be  fitted  to  these  posi¬ 
tions  using  linear  regression  analysis.  This  provides  an 
initial  guess  for  the  phase  function  estimate  Pix, ). 
Linear  regression  analysis  may  also  be  used  with  the 
position  and  magnitude  of  the  fringe  maxima  and 
minima  to  find  initial  guesses  for  the  background  and 
envelope  function  estimates  fil.r,)  and  £(. v,  >,  respec¬ 
tively. 

Many  nonlinear  regression  analysis  algorithms  are 
based  on  linearization  of  the  problem  in  a  local  area. 
Although  an  extensive  evaluation  of  different  algo¬ 
rithms  has  not  been  made,  we  found  that  a  Gauss- 
Newton  ( linearization >  method  from  the  IBM  SSP 
package  failed  to  converge  to  a  valid  solution  in  tests 
with  fringe  data.  It  seems  that  the  irradiance  distri¬ 
bution  may  be  sufficiently  nonlinear  that  algorithms 
based  on  linearization  are  not  useful.  An  optimized 
iterative  search  technique  suggested  by  Powell-  has 
been  used  successfully.  It  is  described  in  some  detail 
in  the  Appendix. 

IV.  Example 

Here  we  present  an  example  which  demonstrates  the 
evaluation  of  a  double-exposure  holographic  interfer- 
ogram  of  a  simply  loaded  aluminum  beam.  A  system 
of  lenses  was  used  to  form  a  full  size  real  image  of  the 
aluminum  beam  from  the  reconstructed  holographic 
image.  A  microcomputer  was  used  to  record  irradiance 
measurements  from  a  photometer  connected  to  a  bO-um 
fiber-optic  probe.  The  probe  scanned  a  I -cm  straight 
path  in  the  plane  of  the  real  image  of  the  beam.  There 
were  197  irradiance  measurements  taken  at  approxi¬ 
mately  uniform  intervals  along  the  path. 

The  particular  data  set  presented  here  is  for  a  region 
which  includes  the  maximum  deflection  of  the  beam. 
The  fringes  are  not  localized  near  the  surface  of  the 
beam.  Fringes  to  the  left  of  the  maximum  deflection 
are  localized  behind  the  surface  of  the  beam.  Fringes 
to  the  right  are  localized  in  front  of  the  beam.  The 
imaging  system  used  a  small  entrance  pupil  so  that  the 
fringes  on  both  sides  of  the  image  are  visible.  This  in¬ 
troduced  a  large  amount  of  speckle  noise  into  the  image. 
Figure  L  shows  the  measured  irradiance  data  tor  this 
image. 

The  positions  of  fringe  centers  tboth  maxima  and 
minima)  were  estimated  visually  from  Fig.  1.  The  po¬ 
sition  irradiance  and  relative  phase  at  the  fringe  centers 
are  listed  in  Table  I.  Linear  regression  analysis  was 
used  to  fit  a  third-order  polynomial  to  describe  the 
phase  of  the  fringes  as  a  function  of  position.  This 
provided  the  initial  guesses  for  the  phase  function: 

Pi  v,  I  a  2. Si*  +  iU.Ojt  -  27.9.x -*  -  2.S2.V  \  ■  •*» 

Figure  2  shows  a  plot  of  this  initial  guess  phase  distri¬ 
bution  together  with  the  fitted  points. 

To  estimate  the  background  and  envelope  functions, 
a  straight  line  was  fitted  through  the  irradiance  the 
left  and  right  fringe  maxima  to  describe  the  amplitude 
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Fig.  1.  Measured  irradiance  vs  position  top  the  fringe  pattern  eval¬ 
uated  in  the  example. 


Table  |.  irradiance  in  Arbitrary  Units  and  Phase  in  Radians  at  Fringe 
C  enters 


Position 

Irradiance 

Phase 

0  01 

9.142 

O.ld 

o.i)2 

0.2S  5 

0.29 

0.17 

9.425 

0.7)0 

0.50 

d 

0.75 

0.19 

9.425 

0.S9 

0.51 

K  2*5 

Of  T 


3  !  _ _ _ _ - _ . 

°:.oo  3/c  :.:c  130  :.40  z.s  a  iso  o.^o  a. so  c.90  voo 

p3SiUan  _  :» 

Fit i.  2.  Initial  guess  tor  the  fringe  pattern  phase  distribution.  A 
third-order  polynomial  was  fitted  to  the  fringe  renters  by  linear  re¬ 
gression  analysis.  The  coefficients  from  this  curve  fit  are  used  as  the 
initial  guess  for  the  nonlinear  regression  analysis. 

of  the  maximum  irradiance.  Another  straight  line  was 
fitted  through  the  irradiance  at  the  three  fringe  minima 
to  describe  the  amplitude  of  the  minimum  irradiance. 
The  average  of  the  polynomials  that  describe  the  max¬ 
imum  and  the  minimum  provide  an  initial  guess  for  the 
background  function: 

Alt,  l  «  0.414  -  <91 

Half  of  the  difference  between  the  maximum  and  the 
minimum  is  the  initial  guess  for  the  envelope  func¬ 
tion: 

£\xt)  -  0 .22$  -0.o7.Ur.  1 101 

With  these  values  as  initial  guesses,  an  iterative 
nonlinear  regression  analysis  algorithm  was  used  to 


select  the  set  of  coefficients  that  yield  the  minimum  sum 
of  errors  squared.  Figure  ‘I  shows  the  irradiance  dis¬ 
tribution  that  resulted  from  this  process  superimposed 
on  the  original  irradiance  data. 

The  phase  function  as  determined  by  the  nonlinear 
regression  analysis  is 

/Yvf  i  =  J.40  +  .Mi.2x  -  IT.'ijr-  +  1.99.1  *.  1 1 1  \ 

Figure  4  shows  this  phase  function  from  the  nonlinear 
regression  analysis  compared  with  the  initial  guess 
phase  function  from  analysis  of  fringe  center  positions. 
The  maximum  difference  between  the  two  estimates  is 
t/:J  rad. 

In  conventional  fringe  recognition,  fringes  are  con¬ 
sidered  to  be  binary— they  are  either  bright  or  dark. 
Using  this  convention,  the  center  fringe  in  this  example, 
which  represents  the  maximum  deflection  of  the  beam, 
is  considered  a  bright  fringe.  Our  method  based  on 
nonlinear  regression  analysis  is  capable  of  recognizing 
that  only  a  fractional  fringe  shift  has  occurred. 

In  this  example  application,  as  well  as  others  we  have 
studied,  the  numerical  values  of  the  sums  of  squared 
errors  indicate  that  nonlinear  regression  significantly 
improves  the  description  of  the  phase  distribution  over 
the  initial  guess.  The  initial  guess  irradiance  descrip¬ 
tion  has  a  sum  of  squared  errors  of  E  *  5.  If  the 

background  and  envelope  functions  Bix, )  and  £<x(  i  are 
optimized  by  regression  analysis  while  the  phase  func¬ 
tion  Pi: r, )  retains  its  initial  guess  values,  the  sum  of 

o 

1 

o 

-N 
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Position  -  cn 

Fig.  X  Irradiance  distribution  Junction  that  resulted  tr> -m  the 
nonlinear  regression  analysis  superimposed  on  the  measured  irra 
diance  data. 


Fig.  4.  Fringe  pattern  phase  function  that  resulted  trom  *he  non¬ 
linear  regression  analysis  compared  with  the  initial  gue**  pha*e 
function  that  wa>  based  on  iringe  center  |x*sition  data 
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squared  errors  is  reduced  to  £  *  2.71.  After  optimi¬ 
zation  of  the  full  irradiance  description  the  sum  of 
squared  errors  is  reduced  to  E  =  2.13. 

V.  Applicability 

Evaluation  of  the  fringe  pattern  described  in  this 
example  and  several  other  l-D  fringe  patterns  has  led 
to  some  understanding  of  when  this  method  will  and  will 
not  contribute  to  evaluation  of  a  fringe  pattern.  This 
method  is  most  useful  when  there  are  only  a  few  fringes 
in  the  image  and  when  there  are  many  irradiance  mea¬ 
surements  per  fringe. 

If  there  are  more  than  a  few  fringes  in  a  fringe  pattern, 
there  is  probably  sufficient  information  in  the  fringe 
position  data  to  allow  accurate  evaluation  of  the  phase 
function.  However,  even  in  a  fringe  pattern  with  many 
fringes  it  may  be  desirable  to  use  this  method  in  local 
sections  of  the  fringe  pattern,  particularly  if  the  data  are 
noisy. 

If  there  are  only  a  few  irradiance  measurements  per 
fringe,  the  information  content  of  all  the  irradiance 
measurements  may  not  be  much  greater  than  the  in¬ 
formation  content  of  the  fringe  position  data.  In  this 
case  regression  analysis  will  not  provide  a  significantly 
more  accurate  estimate  of  the  fringe  phase. 

For  simple  fringe  patterns  like  the  one  in  this  exam¬ 
ple,  with  200  or  fewer  data  points  and  10  or  fewer  coef¬ 
ficients,  an  LSI  1 1/23  microcomputer  has  been  used  to 
perform  the  calculations  for  this  method.  More  com¬ 
plex  problems  have  required  use  of  a  larger  and  faster 
computer. 

In  closing,  we  note  that  a  variety  of  other  techniques 
have  been  suggested  for  determining  phase  distributions 
from  fringe  data.  The  work  of  Brandt  and  Taylor  is  of 
particular  relevance  as  it  deals  with  problems  of  the  type 
discussed  in  the  present  paper  through  an  iterative 
technique  to  fit  spline  functions. u  Sollid0  fitted  2-D 
surface  deformation  patterns  to  Bezier  polynomials.  A 
referee  has  pointed  out  that  several  schemes  based  on 
linear  regression  analysis  have  been  developed  for  in¬ 
terferometric  testing  of  optical  elements.6*8 


Appmdlx:  Numerical  Method 

We  solved  for  all  the  coefficients  in  8(xl L  £{xt  ),  and 
P(x,)  using  a  nonlinear  regression  analysis  algorithm 
proposed  by  Powells  Powells  algorithm  is  an  opti¬ 
mized  version  of  a  type  in  which  one  variable  at  a  time 
is  changed. 

Consider  the  coefficients  of  &{x,  h  £(x,  ),  and  £(x, ) 
to  be  renamed  as  n  coefficients  xr  of  the  sum  of  squares 
of  errors  function  £.  The  object  of  regression  analysis 
is  to  find  the  values  of  the  n  parameters  xxjc  >, .  .  .  ,xn 
so  that  the  value  of  the  sum  of  squares  of  errors 

Eix \.x  * . xn  1  is  a  minimum.  Consider  xx.x  , . xn 

to  be  directions  in  an  n  -D  space.  These  directions  are 
probably  not  the  optimum  ones  in  which  to  search  for 
the  minimum  of  £.  Powell  describes  a  method  and 
rationale  for  selecting  vector  combinations  of  parameter 


directions  that  can  be  expected  to  be  better  for  search¬ 
ing  for  the  minimum  of  £.  He  also  provides  a  criterion 
for  selecting  which  search  direction  to  replace  with  a 
new  search  direction  and  a  criterion  for  deciding 
whether  to  actually  make  the  replacement. 

An  iteration  of  Powells  algorithm  starting  from  point 
po  consists  of  the  following  five  steps: 

Step  I:  A  search  along  n  direction  vectors  X.  For 

r  =  1,2 . rx  calculate  Ar  so  that  £  ( pr_  l  +  K  Xr  )  is  a 

minimum  and  define  pr  *  p r~i  +  ArXr. 

Step  2:  Determine  which  search  gave  the  largest 
reduction  in  £.  Find  the  integer  m .  I  <  m  <  n .  so  that 
(£*(  p^  _  t )  -  E ( pm )}  is  a  maximum,  and  define  A  = 
(£  I  Pm  -  1 )  -  £lP  *>]. 

Step  3:  Calculate  E:\  =  E(2pn  -  p0>  and  define  £ ! 
=  £lpo)  and  E>  =  £(pn). 

Step  4:  Determine  whether  to  replace  the  vector  X^, 
with  a  new  vector.  If  £3  >  £1  and/or  (Ex  *2 £  j  +  £  <) 

•  1  Ei  -  E  >-  A)-  >  -i  £  1  -  E:\Y*/ 2,  use  the  old  directions 

Xj.Xj . Xn,  for  the  next  iteration  and  use  pn  for  the 

next  po  otherwise. 

Step  5:  Calculate  the  direction  for  a  new  vector  X. 
Define  X  =  ipn  -  Po)  Search  along  this  new  vector  for 
the  new  minimum.  Calculate  A  so  that  £(pn  +  AX)  is 
a  minimum.  Replace  the  vector  Xm  with  X^  *  X  for 
use  in  future  iterations  and  let  p<>  =  pn  +  AX  he  the  new 
starting  point  for  the  next  iteration. 

The  search  directions  for  the  first  iteration  are  the 
parameter  directions  themselves. 

Although  the  Powell  algorithm  has  been  found  to  be 
satisfactory  for  this  application,  there  is  no  sense  in 
which  it  is  known  to  be  the  best  method  for  fringe  in¬ 
terpolation. 

The  coefficients  of  £(x, )  and  £(xt )  are  linear  with 
respect  to  the  irradiance  distribution  function  /lx.  ).  It 
should  be  mathematically  much  more  efficient  to  find 
values  for  these  coefficients  with  linear  regression 
analysis.  Only  the  coefficients  of  P(xt )  should  be  found 
with  a  nonlinear  regression  analysis  algorithm  such  as 
that  of  Powell. 


This  work  was  sponsored  by  the  C.S.  Army  Research 
Office. 
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Tomographic  reconstruction  of  strongly  refracting  fields  and 
its  application  to  interferometric  measurement  of 
boundary  layers 


Soyoung  Cha  and  C.  M.  Vest 


An  iterative  algorithm  for  tomographic  reconstruction  of  refractive-index  fields  from  measured  values  of 
path  integrals  along  rays  which  have  been  bent  by  retraction  is  presented.  The  behavior  of  the  algorithm 
is  studied  by  applying  it  to  path  length  data  obtained  by  computer  simulation  of  experiments  in  which  holo¬ 
graphic  or  Mach-Zehnder  interferograms  of  the  field  are  recorded  for  several  different  viewing  directions. 
A  special  form  of  the  algorithm  is  also  used  to  measure  concentration  profiles  in  the  boundary  layer  formed 
at  the  cathode  of  an  electrolytic  cell  containing  ZnCU  The  Appendix  contains  a  discussion  ol  series  expan 
sion  techniques  for  reconstructing  object  fields  from  measured  values  of  line  integrals  through  the  field. 


I.  Introduction 

The  problem  of  reconstructing  a  refractive-index  field 
from  measured  values  of  path  integrals  through  it  is  of 
interest  in  optical  interferometry1  -  and  computer- 
assisted  ultrasonic  tomography.*'*  In  this  paper  we 
address  this  problem  in  the  context  of  measurement  of 
fields  in  which  refraction  causes  significant  bending  of 
the  probing  rays.  Such  strong  refraction  abrogates  the 
direct  applications  of  analytical9  and  algebraic10  re¬ 
construction  techniques  normally  used  for  computer- 
assisted  tomography  because  they  are  based  on  the  as¬ 
sumption  of  straight  probing  rays.  We  present  here  an 
iterative  scheme,  outlined  in  a  previous  Letter,11  for 
reconstructing  strongly  refracting  refractive-index  fields 
from  data  obtained  by  multidirectional  holographic 
interferometry.  The  technique  is  based  on  ray  optics 
and  involves  use  of  an  imaging  system  during  formation 
uf  the  interferogram.  Other  authors  have  recently  re¬ 
ported  iterative  schemes  based  on  ray  optics  and  alge¬ 
braic  reconstruction  techniques  for  applications  to 
computer-assisted  ultrasonic  tomography.4  *  The  ul¬ 
trasonic  techniques  reported  to  date  differ  substantially 
from  the  optical  methods  discussed  herein  because  no 
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imaging  is  used.  In  fact,  some  of  the  most  detailed  work 
dealing  with  strong  ultrasonic  refraction"  involves  the 
use  of  a  transmitter-receiver  pair,  which  scans  laterally 
with  the  two  elements  fixed  relative  to  each  other  as  if 
mounted  on  a  rigid  bar.  McKinnon  and  Bates1-  have 
shown  that  this  data  recording  format  can  result  in  large 
areas  within  the  object  being  inaccessible  to  tomography 
when  refraction  is  significant. 

In  multidirectional  interferometry,  either  a  holo¬ 
graphic-  or  Mach-Zehnder1*  interferometer  is  used  to 
record  interferograms  of  some  object  field  from  a  variety 
of  viewing  directions.  Path  length  differences  are 
measured  by  analyzing  the  fringe  pattern  formed  when 
an  optical  wave  that  has  passed  through  the  test  region 
interferes  with  some  reference  wave.  In  the  case  of 
double-exposure  holographic  interferometry,  this  ref¬ 
erence  wave  is  one  that  passes  through  the  same  test 
region  in  the  absence  of  the  object  field.  In  the  case  of 
Mach-Zehnder  interferometry,  the  reference  wave  is  a 
plane  wave,  which  spatially  bypasses  the  object  field 
before  being  recombined  with  the  object  wave.  The 
following  discussion  applies  to  either  case. 

Figure  1  is  a  schematic  diagram  showing  the  forma¬ 
tion  of  an  interferogram.  An  imaging  system,  repre¬ 
sented  here  by  a  single  thin  lens,  is  used  to  form  a  real 
image  of  a  plane  in  the  object  defined  by  the  coordinate 
r}  .  The  object  refractive-index  field  is  assumed  to  he 
within  a  circular  region.  Two  rays  intersect  and  in¬ 
terfere  at  a  typical  point  P'  in  the  image  plane,  (hie  i> 
the  straight  ray  DEF  through  the  undisturbed  medium 
with  refractive  index  n,>.  This  rav  is  defined  hv  its 
distance  p  from  the  optical  axis.  The  other  ray  .-ffi  is 
curved  due  to  refraction  by  the  refractive-index  field 
n(r,0).  This  ray  is  located  by  projecting  the  emerging 
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Kig.  1.  Formation  *»t  an  interten^nm. 


segment  of  ray  Ah  hack  along  its  tangent  to  the  object 
plane,  thereby  defining  the  apparent  object  point  P.  If 
the  light  entering  the  object  field  is  collimated  and  if  the 
imaging  system  is  ideal,  no  relative  changes  of  optical 
path  length  between  the  two  rays  occur  to  the  left  of 
points  A  and  D  or  to  the  right  of  points  C  and  F.  Point 
C  is  the  intersection  of  the  exiting  refracted  ray  with  a 
circular  arc  centered  at  P  and  passing  through  F.  The 
optical  path  length  difference  of  the  two  rays  which 
interfere  at  point  P'  is  therefore 

ATe-'*’  «  vir.  .  »»/*  f  '»  *  hi  —  / > b.  -  bb  >  •  1  • 

The  first  term  in  Eq.  1 1  Ms  the  integral  of  the  refrac¬ 
tive  index  along  the  curved  path  Ah.  which  is  a  solution 
of  the  ray  equation: 


The  optical  path  length  difference  is  a  nonlinear  integral 
transform  of  nir.o),  which  is  defined  by  Eqs.  \  1 )  and  (2). 
We  term  this  the  path  length  transform  and  denote  it 
by  P: 

A'iqj.W  *  P(N<r..a;rTn|.  *  0 

For  each  viewing  direction,  defined^  by  angle  from 
some  arbitrary  reference  direction,  AT  can  l>e  measured 
by  determining  the  variation  of  fringe  order  number  St 
with  />; 

*  \<\S o '  n 

where  \u  is  the  wavelength  of  light  in  the  undisturbed 
medium  with  refractive  index  n,(. 

If  the  refractive-index  variation  is  sufficiently  small, 
the  curved  ray  Ah  in  Fig.  I  will  essentially  coincide  with 
the  straight  line  EF,  and  Eqs.  <  1)  and  Cl)  will  reduce 
to 

S'ptfijh  3  J*  pjir.-.a  -  n,.\dl .  o’o 

P|n»r...a;.»T,i|  *  R|mr.*;»i  -  n„\  *  A't>»/c*U. 

In  Eqs.  i5>  and  <6)  the  overbar  denotes  quantities  de¬ 
fined  for  straight  paths,  dl  is  a  differential  distance 
along  the  straight  line  EF,  and  P  is  termed  the  line- 
integral  transform  which  is  mathematically  equivalent 
to  the  2-D  Radon  transform  R  of  -  nt}. 


In  this  paper  a  scheme  for  approximate  inversion  of 
the  path  length  transform  P  is  presented.  To  date  no 
analytical  inversion  of  P  is  known  for  the  general 
asymmetric  case,  although  such  an  inversion  is  known 
for  the  case  of  strongly  refracting  radially  symmetric 
object  fields  n(r).u  An  extensive  literature  dealing 
with  inversion  of  the  line  integral  transform  exists  be¬ 
cause  it  forms  the  basis  of  computer-assisted  tomog¬ 
raphy. 

In  Sec.  II  we  present  an  iterative  scheme  for  recon¬ 
struction  of  a  strongly  refracting  refractive-index  field 
from  measured  values  of  its  path  length  transform.  In 
Sec.  Ill  we  present  measurements  of  mass  concentration 
boundary  layers  which  were  obtained  by  applying  the 
iterative  reconstruction  scheme  to  data  obtained  by 
holographic  interferometry.  The  Appendix  contains 
a  discussion  of  the  application  of  series  expansion 
methods  to  inversion  of  the  line  integral  transform: 
these  methods  are  utilized  as  part  of  the  iterative 
scheme. 

II.  Iterative  Reconstruction  of  Strongly  Refracting 
Fields 

Appreciable  errors  may  result  if  data  representing  the 
path  length  transform  of  a  strongly  refracting  field  are 
reconstructed  by  inversion  codes  designed  for  the  or¬ 
dinary  case  of  tomography  in  the  refractionless  limit. 
We  have  used  an  iterative  procedure  for  computational 
inversion  of  the  path  length  transform.  The  method 
is  based  on  successive  estimation  of  the  deviation 
function  defined  by 

=  AT< f,.Ui  -  ATqi.*L 

It  is  assumed  here  that  the  domains  of  definition  of  AT 
and  AT  are  the  same.  In  some  applications  this  is  not 
the  case,  and  a  transform  termed  the  T  transform 
TiAT)  must  he  introduced  to  make  the  domains  >>t 
definition  identical.1 ' 

The  iterative  algorithm  is  as  follows11: 

til  Make  an  initial  estimate  of  the  deviation  function 
[)  i  ffjh.  where  i  * 

<  2)  Calculate  the  corresponding  estimate  of  the  line 
inregraJ  transform,  i.e..  the  refract  ion  less  path  length 
transform 

A’T  * *  A  -  D 

i:0  Approximately  reconstruct  the  field  by  compu¬ 
tational  inverse  line-integral  transformation: 

a  -  n„  *  A'P  t. 

U)  Using  computational  ray  tracing,  calculate  the 
path  length  transform  of  the  estimated  field: 

A‘F  =  P|'j  .a:r»«t|.  1  c 

i  r» i  Calculate  a  new  estimate  of  the  deviation  func¬ 
tion: 

I)  tftjo  *  AT  -  AT  •  1 1  • 

<b)  Return  to  step  i2>  and  continue  the  iterative 
procedure  until  the  change  of  some  measure  of  l)  i/'Tt 
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or  of  the  difference  between  two  successive  recon- 
structed  fields,  is  minimized  or  smaller  than  some  pre¬ 
determined  value. 

To  implement  this  algorithm,  we  expressed  the 
line-integral  transform,  which  must  be  inverted  in  step 
(3),  as  a  Fourier  series  within  the  circular  domain  p/R 

<  1: 

A(t *  £  5-  exp i  121 

^  n 

where  #„,„(/;)  is  the  product  of  an  orthogonal  polynomial 
and  an  envelope  function  whose  specific  forms  are  dis¬ 
cussed  in  the  Appendix.  A  finite  number  of  coefficients 
are  calculated  by  fitting  them  to  the  known  discrete 
values  of  the  line-integral  transform  A *1*1  In 

theory  this  could  be  done  by  the  usual  methods  of 
Fourier  analysis.  In  practice,  redundant  data  were 
recorded  in  an  attempt  to  suppress  amplification  of 
input  errors,  that  is,  a  number  of  coefficients  A„„,  less 
than  the  number  of  given  values  of  ATi/>.  A  )  were  cal¬ 
culated  by  solving  the  overdetermined  system  of  alge¬ 
braic  equations 

l  3  5.  1.  I  ll:li 

in  the  least  squares  sense  using  the  subroutine  DLLSQ 
based  on  an  algorithm  due  to  Golub.17  Once  the  coef¬ 
ficients  A..,,,  have  been  calculated,  the  reconstructed 
distribution  of  change  of  refractive  index,  /ir.o)  =  mr.tf) 
-  n„.  can  be  represented  by  a  Fourier  series  within  the 
circular  domain  r  R  <  I: 

| . ;w...in  •  i/p 

Appropriate  functions  and  ir  \  which  are 

line-integral  transform  pairs  are  discussed  in  the  Ap¬ 
pendix  and  in  more  detail  in  the  thesis  of  Cha.1*’ 

We  chose  the  approach  of  series  expansion  of  the  field 
and  its  transform  because  it  is  convenient  for  accurate 
and  efficient  computational  ray  tracing  and  because  the 
fields  to  which  we  intend  to  apply  the  code  have  rela¬ 
tively  smooth  structures.  The  series  expansion  ap¬ 
proach  used  in  our  investigations  yielded  a  fairly  effi¬ 
cient  reconstruction  code,  and  the  matrix  of  the  equa¬ 
tions  for  is  relatively  sparse.  The  only  nonzero 
coefficients  .-U,„  are  those  for  which 

*>i  *  it. ± l ,± j,±  i . 


/  *  ...  4  1  I 

In  step  (4)  of  the  algorithm,  the  path  length  transform 
of  an  estimated  refractive-index  field  must  be  com¬ 
puted.  Because  the  refractive-index  distribution  is 
expressed  in  analytical  form,  namely,  by  a  series  as  in 
Eq.  (14)  rather  than  by  values  at  discrete  points,  we 
found  it  convenient  to  base  the  ray  tracing  procedure 
on  the  use  of  Hamilton’s  equations1'  rather  than  on 
Snells  law  or  direct  integration  of  the  ray  equation. 
Indeed  a  modified  form  of  Hamilton’s  equations  can  be 
obtained  by  decomposing  the  ray  equation  [Eq.  (2)|  into 
an  equivalent  set  of  first-order  equations.  There  are 


advantages  in  terms  of  accuracy  and  computation  time 
to  such  decomposition.19  In  Cartesian  coordinates,  the 
equations  are 

d\ 

=  tan>.  'Plat 

dx 

<m  dn 

- tarry 

d  >  a\  dx 

—  -  *  'l'ibi 

dx  n 

d  t  A £  )  n 

— - * -  .  .  in<  , 

dx  cos> 

where  A^y  is  the  angle  between  the  tangent  to  the  ray 
and  the  x  axis.  A E  is  the  optical  path  length  of  the  ray. 
which,  when  evaluated  as  the  ray  exits  the  circular  re 
gion  r/R  <  1.  is  the  value  of  the  integral  appearing  in  rhe 
path  length  transform.  Eq.  (I).  The  corresponding 
system  in  polar  coordinates  is 

do  l 

—  *  -  tano.  ■  i"  i« 

dr  r 

d\  1  M"  / 

—  * - I-  [r - +  n 

dr  rn  (  ,«r 

dO£)  n 

d  r  it  »s  ^ 

where  «:>  is  the  polar  angle  from  the  x  axis  to  the  radial 
position  vector  of  a  point  on  the  ray.  and  \  is  the  angle 
between  the  tangent  to  the  ray  at  that  point  and  the 
radial  position  vector.  After  testing  several  numerical 
integration  schemes  with  regard  t.  computational  ef- 
ficiencv  for  this  application,  we  chose  to  use  the  sub¬ 
routine  package  DRIVE-*’  to  integrate  Eqs.  i  lb)  or  i 1 7 ». 
A  variable-order  variable-step-size  differential  equation 
solver  of  the  Adams  type  based  on  Gear's  subroutine 
difsi  B-1  was  used.  The  path  length  transform.*  re 
quired  in  step  (4)  were  computed  in  this  manner  for  -ets 
of  equally  spaced  parallel  rays  entering  the  region  r  R 
<\.  Cubic  spline  interpolation  was  used  to  estimate 
A (PipJh  corresponding  to  other  ray  paths  through  the 
field. 

The  path  length  transform  depends  in  part  on  the 
object  plane  selected  by  the  observer  when  the  image 
of  the  interference  pattern  is  formed.  The  importance 
of  the  selection  of  this  plane  is  evidenced  by  the  fact  that 
if  the  object  field  is  radially  symmetric,  the  deviation 
function  can  he  experimentally  reduced  nearly  to  ,vn* 
by  focusing  on  the  center  plane  of  the  object.’*  A 
similar  effect  occurs  in  boundary  layer-type  fields 
around  objects  of  length  L  if  one  focuses  a  distance  L  I 
from  the  exit  plane.--  If  the  object  field  is  asymmei  ric. 
the  observer  should  select  an  object  plane  in  which  no 
apparent  ray  crossing  occurs  and  in  which  the  tri nge 
spacing  across  the  interferogram  is  fairly  uniform. 

During  operation  of  the  algorithm,  crossing  <>!  com¬ 
puted  rays  through  the  estimated  field  can  occur.  The 
computer  code  should  include  logic  to  recognize  Mich 
occurrences^and  discard  the  corresponding  computed 
values  of  A4>. 

The  behavior  of  the  algorithm  was  studied  by 
applying  it  to  data  generated  by  numerical  simc.  uiun 
of  interferometry  experiments.  Discrete  values  oi  t  he 


15  August  1981  /  Vo!  20.  No.  16  APPLIED  OPTICS  2789 


Error 


Fig.  2.  Reconstruction  of  an  asymmetric  object 
field:  <a)  plot  of  the  absolute  value  of  the  retrac¬ 
tive-index  distribution  described  by  Eq.  <!*);  <b> 
reduction  of  error  by  iterative  reconstruction  of  this 
field. 


path  length  transform  were  generated  by  computational 
ray  tracing  through  test  object  fields  expressed  in  ana¬ 
lytical  form.  The  transform  was  computed  for  several 
different  object  planes,  from  which  one  was  selected 
using  the  criteria  noted  above.  After  each  iteration  the 
current  estimate  of  the  object  field  was  recorded,  and 
the  maximum  error  and  average  error  over  the  field  were 
computed  as  percentages  of  the  maximum  refractive- 
index  change.  In  all  cases,  n0  3  1. 

Figure  2(a)  is  a  plot  of  the  absolute  value  of  an 
asymmetric  object  field  expressed  as 

n  -  *  -0.06(1  -  \r  R)~]\\x  R)  +  l.5|2,  » 18) 

The  simulated  data  consisted  of  thirty-one  equally 
spaced  values  of  A<t>  for  each  of  eight  viewing  directions 
spaced  at  equal  angular  intervals  over  a  180°  range  of 
viewing  directions.  This  field  refracts  quite  strongly, 
bending  some  rays  by  as  much  as  27°.  The  reduction 
of  maximum  and  average  errors  by  the  iterative  proce¬ 
dure  is  shown  in  Fig.  2(b).  In  this  figure,  and  in  those 
that  follow,  iteiation  1  refers  to  reconstruction  by  as¬ 
suming  that  *  A<f>,  i.e.,  refraction  is  neglected. 

Figure  3  shows  the  results  of  reconstructing  a  radially 
symmetric  distribution  which  has  a  single  maximum 
and  decays  gently  to  n0  at  r  *  R: 


OF  iTE*AH0NS 

fb) 

ntr)  -  n«  *  |l  -  0.9  exp(-4(r//?)-||1 

Figure  4  shows  the  results  of  reconstructing  a  radially 
symmetric  field  which  has  two  local  maxima  and  which 
decays  rather  abruptly  to  no  at  the  boundary 

fiir)  -  n0  *  — «Uin[l  —  (r  Rl-|  +  0.025 (cosiOirr  R\  +  l|.  O** 

Although  our  study  of  performance  of  the  algorithm 
was  limited  in  scope,  some  conclusions  can  be  drawn. 
Care  must  be  exercised  in  selecting  parameters  such  as 
the  position  of  the  object  plane  irf  in  Fig.  L).  the  number 
of  terms  in  the  series  used  to  represent  the  object  field 
and  its  transform,  and  the  number  of  rays  traced  to 
calculate  estimates  of  the  path  length  transform.  The 
principal  factor  in  determining  computation  time  is  the 
number  of  terms  in  the  series,  because  the  series  must 
be  evaluated  many  times  during  ray  tracing.  One  way 
of  reducing  the  number  of  terms  is  to  use  an  envelope 
function  whose  form  is  generally  similar  to  the  expected 
form  of  the  field  being  measured.  This  concept  is  de¬ 
scribed  in  the  Appendix,  and  its  application  is  illus¬ 
trated  in  the  next  section.  Although  we  have  demon¬ 
strated  successful  reconstruction  of  fields  with  multiple 
maxima  and  quite  steep  gradients,  there  are  surelv 
limitations  to  the  complexity  of  fields  which  can  be  ac¬ 
curately  reconstructed  using  limited  data  and  com- 
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tractive -index  distribution  described  by  Kq. 

•  a i  field  and  its  reconstruction;  ibi  reduction  rrmr 
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puting  power.  Such  limitations  are  probably  set  by 
excessive  ray  crossing  and  oscillation  of  estimated  ray 
paths  near  closely  spaced  peaks  and  valleys  in  the  re¬ 
fractive-index  distribution. 

III.  Application  to  Measurement  of  Mass  Transfer 
Boundary  Layers 

A  specialized  computer  code  based  on  the  algorithm 
presented  in  Sec.  II  was  written  for  use  in  measuring 
strongly  refracting  boundary  layers.  It  was  applied  to 
interferometric  measurements  of  the  concentration 
boundary  layer  on  the  cathode  of  an  electrolytic  cell 
filled  with  ZnCl?.  Figure  5  is  a  schematic  diagram  of  a 
test  cell  in  which  a  boundary  layer  is  formed  on  the 
surface  of  some  opaque  object.  The  refractive  index 
within  the  boundary  layer  varies  rapidly  and  mono- 
tonicallv  in  the  direction  normal  to  the  object  surface 
but  does  not  vary  appreciably  parallel  to  the  surface. 
The  refractive  index  of  undisturbed  fluid  in  the  ceil  is 
the  refractive  index  of  air  outside  the  cell  is  na ;  and 
the  refractive  index  of  the  window  is  nH . 

We  assume  that  the  interferogram  is  formed  by 
two-exposure  holography.  During  the  first  exposure 
the  fluid  has  a  uniform  refractive  index  flu,  and  during 
the  second  exposure  the  boundary  layer,  whose  refrac¬ 
tive  index  n(.t )  is  to  be  determined,  is  present  adjacent 
to  the  opaque  surface.  The  thickness  of  the  boundary 
layer  is  <v  Optical  rays  which  are  near  the  outer  edge 
of  the  boundary  layer,  x  *  oh ,  are  nearly  undeflected. 
but  those  which  enter  adjacent  to  the  surface,  x  =  0,  are 
strongly  refracted.  Hence,  in  the  interferogram,  the 
object  surface  will  appear  to  be  displaced  as  shown  in 
Fig.  5.  The  amount  of  displacement  depends  on  the 
choice  of  object  plane.  If  the  opaque  surface  is  physi¬ 
cally  present  during  both  holographic  exposures  the 
observer  must  focus  on  the  one  object  plane  in  which 
both  images  of  the  surface  appear  to  coincide:  however, 
apparent  ray  crossing  or  diffraction  patterns  may  make 
the  interferogram  difficult  to  read  near  the  object  sur¬ 
face.  An  alternative  is  to  remove  physically  the  opaque 
surface  during  the  first  holographic  exposure.  In  this 
case,  which  is  optically  equivalent  to  Mach-Zehnder 
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Fig.  V  Refraction  of  an  optical  rav  by  a  boundary  layer  and  um 
-ed ion  window.  Apparent  displacement  of  the  object  Mirfjce  imi 
an  appropriate  distortion  of  the  boundary  laser  re-ult. 

interferometry,  the  observer  is  free  to  choose  an  arbi¬ 
trary  object  plane  in  which  no  apparent  ray  crossing 
occurs  and  in  which  the  fringes  are  more  numerous  and 
more  uniformly  distributed.  We  assume  this  latter  case 
in  the  following  discussion  since  it  is  general. 

The  path  length  transform  is 

P\'u  vi:n..|  =  A-Kim 

i  H  _ 

Ml  V  His  +  M ~  ,  \i  -  ,  „  ' 
t 

+  ft  ,  l  II  “  'I  \/  l  ”  M  ;  (  (  '  J  1  I 

P|ncv km.,|  *  i  =  u\/|mu  i  -  n.  |.  -  j  j  • 

where  x'  is  the  coordinate  measured  from  the  apparent 
location  of  the  opaque  surface.  P  is  defined  over  a  re¬ 
gion  of  width  o(l,  the  apparent  boundary  layer  thickness, 
but  P  is  defined  over  a  region  of  width  To  apply  the 
algorithm,  the  coordinates  must  therefore  be  streu  htd 
by  applying  what  we  term  the  T  transform: 

/‘l-VRr  >|  *  p-Vlu.v  p),  ■  j  :  i 

where  p  -  KlK-  The  line-integral  transform  was  ex¬ 
panded  in  a  polynomial  series  multiplied  by  an  envelope 
function: 
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button  and  its  reconstruction;  *bi  rc*ducti**n  >*t  errors  in  iterative 
reconM  ok  non. 


Prior  to  experimental  work,  several  computer  simu¬ 
lations  were  carried  out.  One  example  is  the  test 
model: 

ft*  1 1  -  •?  *  -n.oiU erto  J.” \  i  =  l.  ’»:>!  1: 

»i.  »  l..VJ:U  n.;  =  I.INNI-JHT4: 

v  .  *  1 1“  mm  <  .  ~  0; 

u  ...  «  jo.on  mm  »Vs  *  1.0  mm. 

Here  erfc  denotes  the  complimentary  error  function. 
Figure  hi  at  shows  the  iterative  improvement  of  the  re¬ 
construction  when  the  plane  is  chosen  so  that  p  =  2.05. 
Similar  results  are  shown  in  Fig.  7  for  the  case  of  p  =  1, 
i.e..  when  the  two  images  of  the  surface  coincide.  In 
both  cases,  ten  simulated  fringe  readings  were  used,  the 
approximating  series  had  five  terms,  and  a  linear  en¬ 
velope  function  etjr)  *  1  —  x,  o*  was  used.  The  simu¬ 
lations  indicate  good  accuracy  and  convergence. 

Tne  test  cell  used  in  the  mass  transfer  experiments 
was  a  Plexiglas  tank  with  flat  quartz  windows.  The 
dimensions  of  the  cell  were  selected  on  the  basis  of  the 
results  of  simulated  experiments.  During  the  experi¬ 
ments,  natural  convection  was  induced  by  differences 
in  concentration  of  ZnCD  across  the  Tmm  wide  gap 
between  the  cathode  and  a  flat  porous  graphite  dia¬ 
phragm.  The  cathode  was  an  accurately  machined 
graphite  bar,  and  the  anode  was  a  zinc  rod  inserted  into 
the  ZnClj  solution  on  the  other  side  of  the  diaphragm. 
Some  physical  properties  and  dimensions  were 


F»ir-  7.  Reconstruction  of  the  same  boundary  layer  as  in  Fiii  ♦?  imaged 
-oihutp  *  1.0;  no  distribution  and  its  reconstruction. 1  lo  reduction 
<>t  errors  by  iterative  reconstruction. 


Fisj.  s.  lntertero^r.un  >•!  .»  portion  >>t  The  <  .ithodu  boundurv  '  c 
U0  current  deii-*it\  >  m.\  cm-:  •  In  current  detain  i-*m\  ,  m 


u.,  =  I.mrj*V7|.  »i,.  =  1  n,  *  l  4 '>71. 

u  *  "»  mm.  .  .  =  t  mm.  <  =>  ti.  to  n  -  • 

Holograms  were  recorded  on  Agfa  Gevaert  lnK75 
plates  using  a  .50-m\V  He-Ne  laser.  The  first  exposure 
was  recorded  with  a  homogeneous  solution  with  con¬ 
cent  ration  Cu  *  1.5  mole/liter  ofZnClj  in  water  at  jo°(\ 
An  electrical  potential  difference  was  then  applied 
across  the  electrodes,  and  the  second  exposure  was  re¬ 
corded  after  the  cell  had  stabilized  for  5  min.  The  m 
terferograms  were  recorded  by  photographing  the  ho¬ 
lographic  virtual  image  through  an  afocal  telescope 
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formed  by  two  828-mm  lenses.  The  interferograms 
were  then  enlarged  bv  a  factor  of  ^  l  Fringe  positions 
in  the  enlarged  photographs  were  measured  with  a  mi¬ 
crometer-driven  translating  microscope.  Interfero¬ 
grams  recorded  with  cathodic  torrent  densities  of  and 
lo  mA  cmJ  are  shown  in  Fig.  >. 

Lsing  the  computer  code  described  above,  the  re¬ 
fractive-index  distribution  in  the  boundary  layer  was 
measured  for  several  different  current  densities. 
Concentration  distributions  were  then  inferred  by  using 
the  relation 

c  —  r  =  *>j.»>ao»  —  i  ) 

which  was  determined  by  retractometrv.  Typical 
boundary  layer  profiles  determined  in  this  manner  are 
shown  in  Fig.  9.  We  were  able  to  measure  boundary 
layer  profiles  for  current  densities  up  to  mA/cm-. 
This  effective  upper  limit  was  set  by  the  loss  of  ability 
to  resolve  individual  fringes  in  the  background  noise 
using  the  available  optical  system.  Other  physical 
factors  such  as  generation  of  additional  ionic  species  at 
high  current  densities  also  limit  the  accuracy  of  con¬ 
centration  measurements. 
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Appendix 

We  utilized  a  series  expansion  method  to  invert  the 
line-integral  transform,  i.e..  the  2-D  Radon  transform 
of  fields  defined  within  a  circular  domain  |r|  <  ft.  For 
simplicity,  let  fir.o)  denote  the  field  to  be  reconstructed 
and  let  gip.'O  denote  its  line-integral  transform.  The 
objective  is  to  expand  gipJh  in  a  Fourier  series  and  then 
invert  that  series  term  by  term  to  obtain  In 

practice  gipJD  is  obtained  by  physical  measurements, 
in  this  case  by  interferometry.  To  speed  convergence 
of  the  series  and  consequently  decrease  computation 


time  we  found  it  convenient  to  introduce  the  concept 
of  an  envelope  function.  eir/R ),  which  is  a  function 
having  roughly  the  same  radial  shape  as  f(r,0)\  similarly 
e(p/R)  has  roughly  the  same  radial  shape  as gip.tt).  It 
is  not  difficult  to  choose  an  appropriate  e(p/R)  by  ex¬ 
amining  the  data.  We  now  show  how  an  appropriate 
series  can  be  chosen  and  inverted  once  eipIR )  has  been 
selected. 

Let  D  be  the  space  of  all  functions  in  2  D  space  w  hich 
have  support,  are  continuous,  and  are  bounded  in  |r| 
<  R.  Furthermore,  let  Dt  tr  be  the  space  spanned  by 
functions  in  D  multiplied  by  e(r/R)%  which  is  assumed 
to  be  an  even  function  which  has  support,  is  continuous, 
and  Lesbesgue  integrable  in  \r\  <  R. 

Tsing  results  set  forth  by  Ludwig-^  and  the  envelope 
function  concept,  it  can  be  show-n  that  the  necessary  and 
sufficient  conditions  for  to  be  the  2-D  Radon 

transform  of  f(r.o)  in  D,  ir  k,  are  that 

( 1 I  -ft.**)  is  a  function  in  D,.Ui  k,.\\ 

12)  gt/i.M  +  rr)  =  gip.ti): 

(2)  J’.T  f-HgipJWp’*  exp i±imH)dpdtt  =  0  if  m  >  k  for 
arbitrary  k  >  0. 

The  line-integral  transform  of  any  function  in  D. ,,  /<, 
can  be  expanded  in  a  series  of  circular  harmonics  and 
orthogonal  polynomials.  If  we  choose  orthogonal 
polynomials  qnUoR )  with  weight  function  *:u>-  ft  >1  I  — 
ft  )-]J  -  and  utilize  conditions  ( l  it  can  be  shown 
that  the  appropriate  series  is  of  the  form 


\i  * 

The  notation  n  =  m:2  indicates  that  n  increases  in  in¬ 
crements  of  2  from  a  lower  limit  m .  i.e.,  only  about  a 
third  of  the  coefficients  in  the  series  are  nonzero. 

The  inverse  transform  can  be  found  on  a  term-by- 
term  basis  using  the  following  relation  among  the  in¬ 
verse  Radon  transform  R~!.  the  1-D  Fourier  tran>torm 
ft,  and  the  Hankel  transform  of  order  nt.H  .  'Ret. 
24): 


R  V1  "O* 

1 1  =  ;  .  t  r  l  ;  r»j  M  | .  \  _ 

where  /...  ir  i  =  / '  H ...  ,F |c„, 

i/>)|l.  Here  the  definition? 

of  the  transforms  are 

-j: 

I'M  *■/.  *  J  T  V  :  i .  \ 

x  » expi  .  -At 

where  Jnix)  is  the  Bessel  function  of  the  tirM  kind 
order  n.  The  following  three  expansions  and  their  in¬ 
verses  are  typical  results  of  this  approach: 

<  D  Let  the  envelope  function  tM piR)  *  1  on  the  disk 
p/R  <  1.  then  the  expansion  of  the  transform  ** />.'** 
is 

^  A  t», |^j  expi  ' \‘>* 

where  P„(p)  is  the  Legendre  polynomial  of  degree  n . 
The  inverse  transform  is 
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f\r,o  l  = 


tK| l  -  (r/fh-|l  - 


•J  ±m  n 

/r \-* 

lr\- 

iy^-1-v.  i^=*i 

/r\-* 

i 

1  —  l- 1 

1  _  V 

-  ' 

1 

\rI 

w 

If  ■!> 

U/  j 

*  A12J 


x 


exptiim^h 


1A6I 


This  expansion  was  first  presented  by  Cha.i:> 

(1)  Let  the  envelope  function  elp/R)  =  (1  —  ip/R)2\l  2 
on  the  disk  p/R  <  1;  then  the  expansion  of  the  transform 
is 


Equation  (AL2)  has  been  presented  in  a  form,  which, 
although  cumbersome,  is  suitable  for  computation  and 
which  explicitly  indicates  that  this  approach  results  in 
the  envelope  function  (1  -  (r/R)-\l  -. 


In' 


exp*  itm*1*. 


( AT ) 


where  Unip)  is  the  Chebyshev  polyn*  mial  of  the  second 
kind  of  degree  n .  The  inverse  transform  is 


where  is  the  Jacobi  polynomial  of  degree  n. 

This  expansion  and  inverse  were  derived  by  Cor- 
mack.-"’ 

(Jl  Let  the  envelope  function  cip)  =  expl  —  a2p2)  and 
let  /  and  a  be  defined  on  the  entire  plane;  then  the  ex¬ 
pansion  of  the  transform  gl pJ*\  is 


=  expt  —  o -/>-*  ^  .A  t„,  t*\p«  iA9* 

where  Hniap)  is  the  Hermite  polynomial  of  degree  n. 
Its  inverse  transform  is 


= —L  exp* -07;')  L  L  -1 1 i  —  1  *  '  1  >-_2"'ar  O' 

V  T  n  —  •}••«  j 

x  L  ’1 : _ .  o:r-‘*  exp*  o*.  <  Aim 

where 


*  a  ,  *  l . 


< n  >..  *  <i*a  +  l  **a  +  2)  .  .  .  *a  +  n  -  2  na  +  n  -  l  >. 


and  L  i r )  is  the  associated  Laguerre  polynomial  of 
degree  n.  This  expansion  and  inverse  were  derived  by 
Maldonado  and  Olsen.-* 

The  computational  inversions  of  line  integral  trans¬ 
forms  performed  while  computing  the  reconstructions 
reported  in  Sec.  II  were  based  on  a  modification  of  Eqs. 
<A5)  and  (A61.  This  modification  resulted  from  in¬ 
troducing  the  additional  constraint  that  g\Rjb  *  0.  i.e., 
that  the  transform  vanishes  at  p  *  ft.  This  results  in 
the  following  expansion  of  the  transform  #(/>,#); 


*  ft.M  * 


V 


V 


exp* 


i  a  i  n 


The  inverse  transform  is 


'll  —  «r  W>J|l  -  ^  ^  .A t„. fxpi ±im<M 


■??  -  'ii  tm  -  n  +■  1 


*  - 1  it— -i  m  +  1 1 _ 2 


> >  •«*  I  J  ii  f  rn  +  1  lt 
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